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ABSTRACT 

A Gibbs-based approach to partitioning hyperspectral imagery into homogeneous regions is investigated for terrain mapping 
applications. The form of Bayesian estimation, Maximum A posteriori (MAP), estimation, is applied through the use of a 
Gibbs distribution defined over a neighborhood system and is implemented as a multi-grid process. Appropriate energy 
fiinctions and neighborhood graph structures are investigated, which model spectral disparities in an image using spectral 
angle and/or Euclidean distance. Experiments are conducted on a HYpICE scene collected over an area adjacent to Fort 
Hood, Texas, that contains a diverse range of terrain features and tfiat is supported with ground truth. Suitable parameter 
ranges are investigated, and the behavior of the algorithm is characterized using individual and combined measures of 
disparity ^yithin the context of a more general framework, one that supports mixed-pixel processing. 

Keywords: Gibbs-based algorithni, Markov Random Field, partitioning, clustering, spectral mixture analysis, hyperspectral 
imagery, 

1. INTRODUCTION 

In prior work,^ an enhanced method of spectral mixture analysis (SMA) was proposed that is appropriate for terrain mapping 
applications, particularly when the imagery involved is characterized by moderate-to-high scene complexity. This type of 
complexity is defined to exist whenever many ftmdamental materials, called "endmembers," exist throughout a scene, or 
whenever some of the fundamental materials in a scetie have spectra that are sinriilar to each other. In such cases, the use of 
one large set of endmembers for performing SMA can be problematic because of the corresponding high dimensionality of 
the SMA model's constrained coordinate systetn, often called its "simplex." A high dinlensional simplex can lead to 
unreliable endmember selection and imprecise estimates of endmember abundance. In order to improve the reliability of the 
model, we proposed to partition a scene into regions that can be treated as homogenous for the purposes of spectral unmixing. 
Once a scene is partitioned into such regions, then we can think of replacing one large set of endmembers by a number of sets 
containing fewer, and more spectrally distinct, endmembers. Each smaller set of endmembers defmes a simplex that is valid 
for the corresponding partitioned regions, and has a lower dimension than the simplex that accpuptSyfor all the endmembers 
in a scene. 

There are many candidate methods for segmenting hyperspectral imagery, both supervised and unsupervised, that could be 
used to partition a scene into regions. However, one desirable attribute, in the proposed setting, is for the regions to exhibit a 
certain amount of spatial consistency. Our approach uses a Bayesian franiework for image analysis that provides a model for 
data with tvyo spatial coordinates and one spectral coordinate (i.e., an image cube). A general mathematical development of 
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the underlying theory, applying the Bayesian framework in the context of image processing, is documented. Using such a 
framework, the properties of a Markov Random Field (MRF) have been applied to develop a Gibbs-based partitioning 
algorithm for gray-scale imagery with two spatial coordinates, and the Gibbs sampler and Metropolis methods have been 
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used as ways to compute maximum a posteriori (MAP) estimates. ' This development was based on Geman and Geman's 
pioi\eering work in applying the MRF framework to the restoration of gray^scale irnagery.^ 

The Gibbs-based partitioning approach that we investigate here processes hyperspectral imagery in a way that imposes spatial 

consistency on the spectral content of sites in each partition. In our previous work,^ the Gibbs-based partitioning algorithm 
was designed with an energy fiinction based on a. spectral angle disparity metric. This algorithm was shown to be sufF!cie;ntly 
sensitive to partition hyperspectral imagery into spatially-consistent regions that are useful for supporting the enhanced 



method of SMA, mentioned above, widiout the clutter often produced by traditional clustering algorithms, such as KMEANS 

and ISODATA.^'^ However, the algorithm was not able to distinguish between certain cultural features (asphalt parking lots 
and certain asphalt-shingled roofs). The sensitivity could be increased to distinguish between these features by. lowering the 
value of a threshold parameter; however, lowering this value would have the undesirable consequence of breaking up natural 
vegetation features into unwanted partitions. 

Herein, two disparity metrics will be investigated - spectral angle and Euclidean distance. A metric shall be defined for each 
of these measures, and an energy fiinction shall be proposed that can respond to the individual measures, as well as respond 
to the combined measures (described in Sections 2.2 and 2.3). We seek to increase the sensitivity of the algorithm without 
introducing undesirable clutter. As an example,.we would like the algorithm to distinguish between cultural features without 
being overly sensitive to uninteresting vegetation changes. 

We investigate tHe.Gibbs partitioning algorithni implemented as a multi-grid process, where the processing typically 
proceeds in stages from course-to-fine grid resolutions and where the beginning (course) grid is initialized randomly. This 
type of initialization is, in some sense, a worst-case scenario. So the algorithm is totally unsupervised. In practice, one could 
initialize the algorithm with the some type of preprocessing scheme - perhaps even another unsupervised or supervised 
classification method. . 



2, SPECTRAIiSPATIAL PARTITIONING AL^ 

2 J Framework for the Partition Algorithm 

A detailed mathematical description of the overall framework for the method discussed here has been formulated,^ where the 
partitioning of hyperspectral imagery is used as a mechanism to condition a spectral niixing process. In the current effort, we 

further investigate the algorithm for estimating the process , called the "partition process," that identifies regions that 
shall be treated as "homogeneous," in the sense that each region contains some culturally similar phenpnienon (not 
necessarily a single material type). 

X^isa discrete labeling Markov Random Field (MRF) process which associates a label with each site. As in the prior 
work, we investigate a multi-grid approach where a is the "grid resolution," which determines the spatial sampling of the 
algorithm at some specific stage of the multi-grid processing. The term "site" refers to a generic element of a lattice. At full- 
resolution a site corresponds to a pixel, but at a coarser scale it corresponds to a square block of pixels. The value of a is 
defined to be the number of pixels per side of each block of pixels. 

We denote S j as the set of sites on an image lattice at full resolution and as the set of sites on a label lattice at the grid 
resolution a: 

At full resolution = 5/ . The label lattice has associated graph structures which describe, for a given element, its 
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neighboring elements. These specified structures allow us to sequentially formulate well-defined probability models for X 
as a Markov Random Field (MRF). IfS isone of S/ or 5 = {s:i ,...,s^}, then an associated neighborhood system 

S = e S} is specified which is a set of lattice graph structures of neighbors, such that S^^g, and 5 6 g). if and only 

ifre^s^ , ' 

The structure of the neighborhood system forms a rather simple pattern. A neighborhood of a site will consist of near, 
intermediate, and far neighbors, and they are all at specific nriultiples of c. Our local neighborhood is either the four closest 



(above, below, left, and right), or eight closest (the surrounding perimeter) sites enclosing the site s. The intermediate and far 
neighbors form a similar pattern as the near neighbors but are at a further distance that is a specific multiple of o. The 
intermediate and far neighbors are used to aid in obtaining a faster global solution and to enable the algorithm to remember 
label assignmentis computed by the algorithm at prior stages. 
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We assurne a common state space where the realizations of are jjCj e F, where F = {1, 2, . . . , I^iabeis } » ^Labels 
the total number of labels assigned to the partitions. The configuration space is then denoted by: * 

af^={x^ =(xf,seS^^^)\xf^{K..,NLMs}} (2) 

is ari MRF with respect to the graph {s]^\ S^^^"^} , where ^^^^^ is the entire neighborhood structure for the 

partition process at resolution a , and ^^^^ is the specific neighborhood at site s for the partition process at resolution O . 

Within this setting, our goal is reduced to finding the Value of X which maximizes ?r(X \ G) , for the observed G, where 
. G i$ the hyperspectral image cube. 

We use the Gibbs equivalence theorem,^ and with an appropriately defined graph , we can develop a Gibbs- 

based algorithm. The advantage of the Gibbs form is that we can compute the MAP estimate, an estimate that maximizes the 
probability Pr(x^ | g) by iteratively sampling from the local Gibbs distribution pertaining to each site sg S^p ^\ This local 
distribution is defined as 

I -yUsixf.g) ' . r1^ 

= V-^ • " (3) , 
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where ^si^s is the energy interaction of site 56 S^p^ with the neighborhood S/*^^^ . We 

use the symbol g , with the underscore, to emphasize the spatial and spectral nature of the realizations of G, which is in 

contrast to the spatial-only realizations of the simpler G.^ Fortunately, because the approach we use to compiute the MAP 
estimate only requires knowledge of the ratios of (3) for current and proposed configurations, the burdensomb computation of 
does not have to be performed. 

We now discuss the specifics of the implementation, below, in Sections 2.2 to 2.4. We define a set of disparity metrics, a 
number of specific energy fimctions of the local Gibbs distribution defined by (3), and then discuss a method of computing 
the corresponding MAP estimates. 

2.2 Disparity Metrics 

The energy functions that will be introduced in Section 2.3 depend on distance measures, called "disparity metrics," which 
indicate the dissimilarity between a site s and a neighbor r. We shall investigate disparity metrics that measure distance in 
terms of the spectral angle metric and the Euclidean distance. 

For a piair of sites at the site locations r and s, consider a set Z)^ ^ = { , D^^] } , where the elements of Dj- J are the 

spectral angle and the Euclidean distance Dpj , between the sites s and r. The individual distance measures 
defined as follows: ' 



2.3 Energy Functions 

Energy functions shall now be defined that respond to disparities as measured by the spectral angle metric, the Euclidean 
distance and perhaps others. As mentioned earlier, we seek to increase the sensitivity of the algorithm without introducing . 
undesirable clutter. 

For each individual disparity metric in (4) ^ (5), an energy function 

wherethedelta fonction ^(jc^.jc^)^ 1 for jc^ = and ^(jc^^ . 

The function f^{Dj!\) maps to the set {-I, I), where = - 1 indicates a npn-significant distance between r 

and s, and a value // ) ^ 1 indicates a significant distance. If < f^^j^^^^^ > then /^^^r'l j = -1 » ai^d otherwise 

, where the threshold ^^^.^g^ determines the sensitivity of the algorithm to the Z)^[^ distance measure. The 
effect of this energy function is to accjumulate the value of /^i^^jl j whenever the label of a site s is equal to its neighbor at r. 

For i=J,2, the value of the threshold ^^f^-^^fj is determined through observation; specifically, K^j^^^^ is determined by 

observing the differences in spectral angle for samples of the same material verses the differences in spectral angle for 

samples of different materials, and ^^^^^ is determined by observing the differences in Euclidean distance for samples of 

the same material verses the differences in Euclidean distance for samples of different materials. If the imagery is calibrated 
to reflectance, then the value attained should be fairly universal, presuming a representative set of materials has been . . 
analyzed. If the imagery has not been calibrated, then this value will depend on the range and scaling of the uncalibrated 
data, 

Iri addition to the energy ftmctions (6) that respond to a sirigle disparity metric, we can construct an energy function that 
responds to a pair of disparity metrics: 

Vs{4,g)^ J,^i4\4y^^AMo^^^^^ (7) 

2.4 Calculating the MAP estimate ' 

A simulated annealing technique is used to obtain the MAP estimate via the Metropolis algorithm, ' ' with an eventual 
fast freeze at low temperature. The computations involve a temperature at iteration k that defines an annealing schedule. 

We use a commonly applied annealing schedule,^'** namely, . 



c 



(8) 
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C is a parameter that corresponds to the maximum depth of a *Vell" in the energy function (6) or (7), The setting for C 
affects the rate of convergence. For example, an upper bound on this value for some disparity / can be computed as the value 

of (6), for which the maximum distance fj Wr}s ) » site-neighbor pairs are all different, but the labels are assigned the 

same value. When running the multi>grid algorithm, computing time can be reduced considerably by using a large C value at 
' the coarsest grid, and otherwise using smaller C values at the fmer grids. Practical values of C are discussed in Sections 3 
and 4. v ■ ^ 

3. DESCRIPTION OF EXPERIME 

IhQata " — ■ ' [ ■ . ■ /. ' 

A HYDICE scene is used iri the experiment that was extracted from data collected over a geographic area near Fort Hood, 
Texas at a 10,000 feet altitude, on 16 May 1996. The scene consists of 300 samples by 600 lines of pixel vector data with a 
spatial resolution of approximately 3 meters. Complete HYDICE image cubes have 210 bands in the approximate 0.4 ^2.5 
\im spectral range. However, because of the severe degradations that occur in atmospheric windows of the spectruni and 
sensor degradations occurring in certain bands, a maximum of 1 17 bands were used to compute the.disparities within the 
energy functions. These 1 17 bands of data were calibrated to reflectance by the Empirical Line Method based on the known 
spectra of the calibration panels that were located in a region outside this particular scene. Ground truth, describing certain 
ground features in the scene, is based on data collected during a site Visit to the area during April 2000. , Figure I shows a 
ground-truth overlay superimposed over Band 49 extracted from the HVDICE image cube. Descriptive information . . 
pertaining to the ground truth is listed in Table I. 
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Figure t. HYDICE Scene (Baad 49) with Ground-Truth Overlay 



Because of the time difference between scene acquisition and the site visit, not all portions of the scene could be documented 
with certainty. In addition to the ground-feature data, spectral data of a diverse range of materials, measured by field 
spectrometers and collected independently of this experiment at other geographic locations, were used to determine an 
appropriate quantization threshold for the spectral angle and Euclidean distance disparity metrics, as discussed below. 
Ground-measured spectral data within the study site were not available. 

3.2, Validation of the Algorithm 

An experiment was conducted using computer code written in C++ that implements our algorithm in the manner described in 
Section 2. The energy functions in (6) and (7) were implemented with the disparity metrics defitied in (4)-(5). In order to . 

determine appropriate values of the thresholds K^^^^^^ and ^^^^^^^ the spectral angles and Euclidean distance measures 
between all the pairwise combinations of the spectra in the reference spectral library, just mentioned above, were computed. 

Tfl|ble 1. Description of the ground truth (GT) information, as displayed in Figure 1. 



Gt Label 


DescriDtion of the Area 


Ai 


Asphalt parking lot, newlyjpaved. 


A2-A3 


Asphalt parking lot A3a and A3b opposite sides of apartment complex, different ages. 


A4-A5 


Asphalt road 


Bl 


Rooftop of a department store, light-gray asphalt shingles, skylights, yellow gas pipes. 


B2a,b 


School, rooftop tnostly coated bubbly light-gray material (i$2b), top edge is gravel (B2a). 


CI 


Concrete pavement to loading dock, some tire marks, 


' ; C2a,b 


Concrete parking lot surrounding small bank building. C2b is brighter than C2a. 


Dl 


Trees: Black Willow, Texas Sugar Berry, Dogwood, Texas oak, Elbon bush. Green Bryer. 


P3 


Trees ; . . 


D4 


Trees: Texas Oak, Burmelia, and Texas Hacfcberry. 


P5 


Trees: Red i3ud, Texas Oak, and Red Bud. , - - 


D6a,b,c 


Trees - Deciduous holly, Bois D'Arc, Osage Orange, Texas Oak, and Bow Wood 


D7 


Scattered junipers throughout the grassy areas. 


El 


Church and school with metal rooftop. 


E2-E6 


Corrugated steel buildings. 


G1-G6 


Tall wild grass. G5 and G6 near drainage and healthier. 


HI 


Healthy maintained grass near bank building. 


Jl 


Private residence, stand-alone. 


J2 


Apartment buildings with asphalt-shingled rooftops. 


J3 


Church with dark brown asphalt-shingled rooftop. 


LI 


Shade of Apartment building. 


Ml 


Tennis court, concrete with blue/gray coating, and lots of dirt.' 


Nl 


Playground area with exposed soil. 


01 


Encircled area is residential with moderately-priced homes and concrete driveways. 


Rl 


Light linear pattern is exposed bedrock. 


R2-R3 


Exposed soil. 


VI 


Asphalt road intersection with exposed soil along the shoulders. 


Wl 


Playground - grass with exposed soil, gravel, and some asphalt. 



The experiment consisted of two trials. Trial A used the spectral angle disparity metric in the energy fimction and Trial B 
used a combined spectral angle and Euclidean distance metric in the energy fimction. The behavior of the algorithm is 
affected by various parameters; namely, o (the grid resolution), C (controls the annealing schedule through the temperature 

)' total (^*"^^ number of iterations through the image), and ^ labels (*^^ number of partition labels). Table 2 lists the 
values of these parameters. As indicated in Table 2, the trials proceeded in a multi-grid sequence with grid resolutipiis 



O = 1 6, 8, 4, 2, 1 . The unsupervised ISODATA clustering method is used as a reference method, 
of ISODATA on the data for 6 classes. 



Figure 2 shows the results 



Table 2. Control parameters for Trials A and B. The sequence of numbers in the "Dist of neighbors" 
column corresponds to the distance from the center pixel (the "site*') to the perimeter, for the surrounding near 
neighbors, the intermediate neighbors, and the far neighbors. The sequence of nunibers in the "# of neighbors'* 
colunin corresponds to the number of near neighbors, intermediate neighbors, and far neighbors in the 
neighborhood system. 



a 


C 


# iterations 


Dist. of neighbors 


# of neighbors 


16 


9.0 


50,000 . 


16,64,128 


8,8,8 


8 • 


'3.0 


3,000 


8,32,64 


8,8,8 


4 


3.0 


3,000 


4,16,32 


8,8,8 


2 


2.0" 


700 


2,8,16 • 


8,8,8 


1 


1.0 


350 


1.4, 8 ■ ' 


8, 8,8 



During the trials, a range C = (1 .0-9.0) was investigated, which is below the upper bound discussed in Section 2. In general, 
the value of C was set to its highest value at the, initial stage of a run (corresponding to high temperature) and it was set to 
lower values for the subsequent stages. In addition to the practical cpmputatipnal advantage, a lower C value allows the 
algorithm to "remember" the results at previous stages. The distances of the intemiediate and far neighbors shown in Table 2 
corresponds to fixed multiples of the grid resolution, ^cr and <9cr, respectively. For all the trials, the algorithm was initialized 
randomly^ and then proceeded froni a grid resolution of (1^/6 to 

4. RESULTS 

Figure 2 and Table 3a shows the results of the ISODATA clustering method applied to the image data for 6 classes. 
Although ISODATA provides a fairly good partitioning of the scene, it is overly sensitive to minor vegetative differences in 
the forested areas, denoted by "D1-D6" in the ground-truth map (Figure 1) as well as by shade. In particular, Table 3a shows 
a very important problem, where another label associated with Al (asphalt) is prevalent throughout subregions of P1-D6 that 
can be attributed to shade. Specifically, 27% (1015 out of 3764) of the D1-D6 pixels are labeled with the same cluster 
associated with asphalt. Also 15.7% (559 out of 3559) of the GKG4 pixels are labeled with the same cluster as associated 
with trees. Furthermore, the clustering does not provide a smooth and spatially consistent mapping of regions, and exhibits a 
significant amount of speckling. 

Based on observations of spectral angle distances between the various pairs of materials of our previous work,^ we set 
^^thresh^ 1 1 .0, which remained constant throughout the partition trials. Based on observations of Euclidean distances 

between the various pairs of materials, we set 1 00.0, which also remained constant throughout the partition trials. 

For both metrics and most of the pair-wise signature combinations, there is an obvious gap between similar materials and 
different materials, except for combinations involving a vegetation signature. Unfortunately, vegetation has a fairly large 
variance, which makes the threshold boundary a bit ftizzy. During these trials, we observed that the exact value of these 
thresholds did not seem to be very critical. 

Figure 3 shows the results of Trial A testing the spectral angle metric in the energy function. On the positive side, the 
algorithm with this energy function was successful in characterizing the structure of many of the natural and cultural features 
in the scene without much of the undesirable clutter prevalent in the ISODATA scene shown in Figure 2. For example, the 
algorithm was not overly sensitive to shade differences such as in the forested areas, denoted by "DI-D6" in the ground-truth 
map. Unfortunately, the algorithm was also insensitive to other differences that we would like to discriminate. For example, 
there is missing structure within the road intersection denoted by "VI,". and the rooftops of the apartment buildings (asphalt 
shingles) denoted by "J2" are confused with the surrouiiding asphalt parking lot "A3." 



Table 3. Confusion matrices of results. Labels in the 18 rows are the GT labels, and labels in the .6 columns are the partitioned 
labels. As an example of how to read the tables: In Table 3a, the region of ground truth extracted for C 1 , consisted of 1 20 pixels. Of the 
1 20 pixels, 1 3 pixels were assigned to ISGDATA cluster label lSO-04, and 107 pixels were assigned to. ISODATA cluster label ISO-05. 

Table 3a. Results for the ISODATA algorithm. 



Labels 


ISO-01 


ISO-02 


ISO-03 


ISO-04' 


ISO-05 


I$O-06 


Al 


0 


0 


518 


0 


0 


0 


Bl 


0 


0 


0 


.278 


0 


0 


Ct 


0 . 


0. 


,0 


13 


107 


0 


C2a 


.0 


0 


0 


0 


13 


0 


C2b 


0 


Q 


.0 


0 


11 


0 


B2a 


^ . . 0 


0 


0 


2 


15 


0 


B2b 


Q 


0 


0 


12 


no 


0 


Ml 


0 


0 


0 


43 


0 


0 


J3 


0 


Q 


33 


44 


0 


0 


A2 


0 ' 


4 


0 


51 


0 


0 


A3a 


0 


0 


0 


16 


0 


0 


A3b 


0 


0 


0 


12 


0- 


0 


J2 


Q 


0 


84 


0 


0 


0 


Pl^ 


' . 2276 


36 


1015 


. 0 


0 


1 


G1-G4 


559 


2989 


3 


0 




8 


HI. 


0 


0 


0 


0 


■ 0 


33 


95 


;lOl 


' Q 


0 


0 


0 


71 


G6 


6' 


IP 


0 


,0 


0 


176 


3b. Results for the Glbbs-based algorithm using the combined disparity (CD) metrics. 


Labels 


CD-05 


CD-01 


CD-02 


CD-04 


CD-03 


CD-06 


Al 


0 


0 


518 


0 


0 


0 


Bl 




0 


0 


278. 


0 


0 


CI 


0 


0 


0 


0 


120 


0 


C2a 


0 


0 


0 


7 


6 


0 


C2b 


0 


0 


0 


11 


0 


0 


B2a 


0 


0 


0 


17 


0 


0 


B2b 


Q 


0 


0 


122 


0 


0 


Ml 


41 . 


0 


0 


0 


2. 


0 


J3 


0 


0 


77 


0 


' 0 


0 


A2 


0 


0 


55 


0, 


0. 


0 


A3a 


0 


0 


0 


0 


,16 


0 


A3b 


0 


0 


0 


0 


12 


0 


J2 


0 


0 


84 


0 


0 


0 


Dl-6 


•3527 


66 


26 


62 


48 


35 


Gl G4 


1 


3543 


0 


I 


0 


14 


HI 


- 33" 


0 


p 


0 


0 


0 


G5 


0 


0 


0 


0 


0 


172 


G6 


0 


185 


0 


0 


1 


0 
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Figure 4 and Table 3b shows the results of Trial B using both the spectral angle and Euclidean distance metrics in the energy 
function. A comparison of Figures 3 and 4 shows a clear advantage to using the two metrics in the energy function. Many of 
the cultural structures missed in Trial A are detected in this trial without introdijcing undesirable vegetation clutter. For 
example, the structure in "Vl" is now evident and so are the rooftops in "J2." In addition, Table 3b shows that the regions 
D1-D6 are labeled quite homogeneously and are not adversely affected by shade, as was the case with ISODATA. 
Specifically, 0.7% (26 out of 3764) of the D 1-D6 pixels are labeled with the same cluster associated with asphalt. Also 0.4% 
( 1 6 out of 3559) of the G I -G4 pixels are labeled with the some other cluster (trees or otherwise). 

The above trials corresponding to Figures 3-4 were performed using identical control parameters, as listed in Table 2, in 
order to provide some consistency for comparing the different methods. The results indicate that an energy fiinction using 
the combined spectral angle and EucUdean distance metrics is the preferred method. 

Overall, the Gibbs-based algorithm generates partitions that accurately represent structured regions in the image. For most 
of the bigger regions, pixels with the same label represent the same type of phenomenon on the ground, globally across the 
image, and regions with different labels correspond to different phenomenon, However, a true globally consistei\t labeling of 
the partitions vvas sonietimes problematic. In particular, labeling problems can be seen in some of the smaller regions. 
Although the contiguous pixels in smaller partitions represent the same ph^^ some small partitions, 

separated by a certain amount of distance and given the same names, actually represent different phenomenon. Conversely, 
some small partitions, separated by a certain amount of distance and given different names, actually represent similar; 
phenomenon. 



5. CONCLUSIONS 

Two energy functions were tested involving the use of the spectral angle metric, and a combined spectral angle metric and 
Euclidean distance, respectively. The results of our experiment indicate that the energy function using a combination of the' 
spectral angle and EucHdean distance metrics produced the best resiiits. ^ 

Overall, the Gibbs-based algorithm generates partitions that accurately represent structured regions in the image. A globally 
consistent labeling of the partitions was sometimes problematic, particularly, for some of the smaller regions. This global 
labeling problem can be significantly reduced by using a non-random initialization scheme, or perhaps by performing more 
than one pass through the multi-grid sequence. 

The Gibbs-based algorithm is envisioned to work within an overall framework involving other processes. In particular, the 
partitioning of hyperspectral imagery is viewed as a mechanism to condition a spectral mixing process. Under this 
fi*amework, the ultimate identification of materials (and the corresponding abundances) in a scene is a consequence of the 
spectral mixing process, so that a precise global labeling is not required at the partitioning stage. 
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A Spectral Mixture Process Conditioned by 
Gibbs-Based Partitioning 
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Abstract — An enhanced method of spectral mixture analysis is 
investigated for hyperspectral imagery of moderate-to-high scene 
complexity, where either a large set of fundamental materials may 
exist throughout, or where some of the fundamental members have 
spectra that are similar to each other. For a complex scene, the use 
of one large set of fundamental materials as the set of "endmem- 
bers'' for performing spectral unmixing can cause unreliable esti- 
mates of material compositions at sites within the scene. In such 
cases, partitioning this large set of endmembers into a number of 
smaller sets is appropriate, where the smaller sets are associated 
with certain regions in a scene. Herein, a Gibbs-based algorithm is 
developed to partition hyperspectral imagery into regions of simi- 
larity. This partitioning algorithm provides an estimator of an un- 
derlying and unobserved process called a ^'partition process" that 
coexists with other underlying (aiid unobserved) processes, one of 
. which is called a "spectral mixing process." The algorithm exploits 
the properties of a Markov random field (MRF) and the associated 
GIbbs equivalence theorem, using a suitably defined graph struc- 
ture and a Gibbs distribution to model the partition process. Con- 
sequently, spatial consistency is imposed on the spectral content 
of sites in each partition. The enhanced spectral mixing process is 
then computed as a linear mixture model that is conditioned on the 
partition process. Experiments are performed using scenes of HY- 
DICE imagery to validate the algorithm, where spectral mixture 
analysis is performed with and without conditioning on the parti- 
tioning process. 

Index Terms — Bayesian, hyperspectral, multigrid Gibbs distri- 
bution, partition process, spectral mixture analysis. 

L INTRODUCTION 

B EFORE introducing our approach, we first give a broad 
overview of the diverse array of detection/classiflca- 
tion/identification techniques which have typically been 
applied to hyperspectral imagery. They are inherently either 
full^pixel or mixed pixel techniques, where each pixel vector 
in an image records the spectral inforniation. The underlying 
assumption goveming full-pixel techniques is that eiach pixel 
vector measures the response of predominantly one underlying 
material/signal at each site in a scene. In contrast, the under- 
lying assumption goveming mixed-pixel techniques is that 
each pixel vector measures the response of multiple underlying 
materials/signals at each site in a scene. Unfortunately, an 
image can often be a combination of the two situations, where 
many sites in a scene are pure materials, but many others are 

Manuscript received November 8, 2000; revised January 30, 2000, 

R. S. Rand is with the U.S. Department of the Army Engineer Re- 
search and Development Center (ERDC), U.S. Corps of Engineers, 
Topographic Engineering Center, Alexandria, VA 22135 USA (e-mail: 
robert.s, rand@erdc.usace.army.mil). 

p. M. Keenan is with tfie Department of Statistics, University of Virginia, 
Charlottesville, VA 22903 USA (e-mail: dmk76@vii^inia.edu). 

Publisher Item Identifier S 0196-2892(01)05503-6. 



mixtures of materials. For perspective, we briefly discuss some 
of the lull-pixel and mixed-pixel techniques, and then introduce 
the enhanced spectral mixture analysis concept. 

By construction, the modeling and algorithms investigated 
here are intended for mapping applications, rather than targeting 
applications. The approach encourages spatial integrity of any 
assigned class labels, so, consequently, the detection of single- 
pixel and subpixel anomalies are discouraged, unless there is 
strong spectral evidence to the contrary. Also, of possible in- 
terest is the following: although the current effort focuses on the 
use of a Gibbs-based partitioning algorithm for enhanced spec- 
tral mixture analysis, this partitioning algorithm can also be used 
with any of the full-pixel methods mentioned, to add spatial in- 
tegrity to an otherwise spectral-only method, as discussed later. 

A. Full-Pixel Techniques 

The simplest and most straightforward full-pixel technique is 
the method of spectral niatching. Spectra of interest in an image 
are matched to training spectra that are either obtained from 
a reference library of spectra or from the image itself Simple 
metrics for determining the degree of match are such measures 
as Euclidean distance, derivative difference, and spectral angle 
[I], [2]. Despite its simplicity, spectral matching can be very ef- 
fective as long as the training data are properly calibrated with 
respect to data of interest and provided that the full-pixel sce- 
nario is appropriate. Unfortunately, a certain number of the pixel 
vectors in a scene will likely measure the spectral response of a 
mixture of materials. If the number of the mixed pixels is signif- 
icant with respect to the size of the scene then spectral matching 
should not be used, at least not as an isolated technique. Also, 
class label assignments generated by spectral matching algo- 
rithms are not affected by spatial neighborhoods, which may 
be a positive consequence for targeting applications. However, 
for mapping applications, consistency of class labels in local- 
ized spatial neighborhoods, termed as "spatial localization," is 
important. 

Other full-pixel methods that have been applied to hyper- 
spectral imagery are the various supervised and unsupervised 
scene segmentation techniques, which are based on statistical 
and pattern recognition methods that were previously used for 
lower-dimensional multispectral processing. As with spectral 
matching, the required training is accomplished with data from 
either spectral libraries or imagery, again with the same caveats. 
The techniques include statistical linear discrimination (e.g.. 
Fisher's linear discriminant), quadratic multivariate classifiers 
(e.g., Mahalanobis and Bayesian maximum likelihood [ML] 
classifiers) [3]-[5], and neural networks [6]. The quadratic 
methods require low^dimensional pixel vectors, and hence. 
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these techniques are typically preceded by a data reduction 
operation to reduce the number of bands. Such neural networks 
as the multilayer feedforward neural network (MLFN) can 
be constructed to model quadratic and higher-order decision 
surfaces without the need to reduce the number of bands, and 
assessments of such methods have been performed [7]-[lO]. 
Studies comparing neural networks to traditional quadratic 
classifiers have also been performed on multispectral imagery, 
with a general consensus favoring the neural methods [11], 
[12]. An advantage of certain neural techniques such as the 
MLFN, is that they can be trained to identify materials that have 
been perturbed by a limited amount of mixing. The methods 
usually do not include any type of spatial localization in the 
decision making process. 

The - most commonly applied unsupervised algorithms 
for clustering multispectral and hyperspectraP imagery are 
KMEANS and ISODATA [4], [5], where the metric used in 
determining the cluster membership is the Euclidean distance. 
Unfortunately, this metric is not particularly responsive to the 
fine structure, or shapes, in high-resolution spectra, and is 
often overly sensitive to intensity differences. Also, these two 
methods dp not include any type of spatial localization in the 
clustering operation. 

B. Mixed-Pixel Techniques 

Spectral mixture analysis (SMA) techniques have overcome 
some of the weaknesses of full-pixel approaches, by using 
linear statistical modeling and/or signal processing techniques 
[ 1 3]-[l 7]. The basic equation is 

x,^Hps + n, (1) 

in matrix notation, where, at site s in the scene, is the pb^ 
served reflected energy, is the modeling parameter vector that 
will be a associated with the mixture proportions, and r/^ is the 
random variable for model error. The columns of matrix H con- 
tain the spectra of the pure materials, called "endmembers." A 
fundamental aspect of basic SMA is that the matrix is pre- 
sumed known and fixed (nonrandom). The assumption of a fixed 
H is problematic, because for most real-world materials, there 
are no single fixed spectral signatures to populate H that can 
precisely represent the pure materials. Further aspects of SMA 
are discussed in Section HI-B. 

Modifications to the basic SMA method often involve parti- 
tioning the matrix H in (1) into desired and undesired signa- 
tures and then computing subspace projections that are orthog- 
onal and/or oblique to the undesired signature and noise compo- 
nents. For example in [18], the concept of orthogonal subspace 
projection (OSP) is applied to hyperspectral imagery that sup- 
presses undesirable signatures and detects a specific signature 
of interest. In this method, H = [D, Lf\, where D is the known 
spectra for a target of interest and is a matrix of undesired 
(but known signatures). In [19], a modification of this technique 
is presented to detect the presence of a known signature in the 
presence of an unknown and undesired background, where the 
columns of U are unknown. Ho>yever, D has to be a minor coni- 
ponent of the scene [19]. These modifications to the basic SMA 



and other similar ones [20] are best suited for targeting applica- 
tions. 

Unless we have accurate ground-truth information about the 
materials in a scene, the task of determining the number of end- 
members in a scene and subsequently identifying these end- 
members is not trivial. Some progress has been made in cases 
where ground truth is not available, but reliable solutions re- 
main allusive, particularly for scenes of moderate-to^high com- 
plexity. Usually, the first step is to "noise whiten" the data and 
then assess its dimensionality. A popular technique is the min^ 
imum noise fraction (MNF) transform [2 1 ]. The MNF transform 
requires estimating a noise covariance matrix J^^, which is typ- 
ically done using the minimum/maximum autocorrelation fac- 
tors (MAF) method developed in [22]. A noise adjusted prin- 
cipal component (NAPC) method was later shown to be math- 
ematically equivalent to the MNF but easier to compute [23]. 
An interactive method using MNF transformed data to identify 
the endmembers has been developed, based on convex geometry 
■ concepts [24], [25] and the minimum volume transform [26]. 
An automated method has been proposed in [27], that of ap- 
plying the minimal descriptive length (MDL) criterion to noise 
whitened data for determining the number of endmembers, fol- . 
lowed by an orthogonal subspace projection method for identi- 
fying them. More recently, methods.of independent component 
analysis (ICA) have been investigated [28], [29]. Unfortunately, 
the reliability of these methods hinges on an ability to accurately 
estimate Difficulties in obtaining reliable estimates of 
are discussed within a targeting framework in [30], in the con- 
text of comparing MNF and PCA transform methods. 

C. Proposed Approach ? 

The present investigation considers hyperspectral image 
analysis from a broader perspective than the individual 
methods mentioned above. Our approach will use the paradigm 
of Bayesian image analysis. A treatment of the various aspects 
of Bayesian image analysis is presented in [31], such as 
the Bayesian paradigm, the Gibbs sampler, and methods of 
estimation. Additional mathematical development of the un- 
derlying theory, applying the Bayesian framework to Random 
Processes in the context of image processing, is presented in 
[32]. Bayesian approaches using Markov Random Field (MRF) 
models have provided algorithms for producing unsupervised 
partition maps and detecting boundaries on single-band im- 
agery [33]. Algorithms by Geman et dl [34] use the Gibbs 
sampler and Metropolis methods to compute maximum a 
posteriori (MAP) estimates. Recently, a Gaussian MRF ap- 
proach that considers the spatial-spectral correlation structure 
of multispectral imagery using the Iterative Conditional Mode 
(ICM) has been developed [35]. Although multiband imagery 
has been treated to some extent in the Bayesian community, the 
significantly different character of hyperspectral imagery has 
not really been addressed. 

In pur experiments, we investigate the Gibbs partitioning al- 
gorithm as an unsupervised niethod, where there is random ini- 
tialization (in some sense, a worst-case scenario). In practice, 
one can consider potential preprocessing methods for initializa- 
tion. One could initialize the proposed algorithm with one of the 
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full-pixel techniques mentioned above, in which case the algo- 
rithm can effectively become supervised. Although this study 
limits itself to random initialization, we see these alternative 
iriitialization schemes as an important topic for future investi- 
gation. In this manner, our approach of incorporating spatial lo- 
calization can potentially improve the results of a number of 
spectral-only techniques. 

An enhanced method of spectral mixture analysis is devel- 
oped that should be appropriate for scenes of moderate-to-high 
scene complexity. This type of complexity is defined to exist 
whenever many endmembers exist throughout a scene, or when- 
ever some of the fundamental members have spectra that are 
similar to each other. For a complex scene, the use of one large 
set of endmenibers for performing SMA can cause unreliable es- 
timates of material compositions at sites within the scene. Our 
approach is to partition such a large set of endmembers into a 
number of sets containing fewer, and more spectrally distinct, 
endmembers, where the smaller sets are associated with certain 
regions in a scene. Recalling the above mentioned fundamental 
(and historical) difficulties in automating endmember selection, 
during the experiments we will be extracting our endmember 
spectra manually, and assigning them to sets based on known 
ground-truth information, in order to control sources of error. 

The remainder of the paper is organized as follows. A math- 
ematical description of the Bayesian model is presented, fol- 
lowed by a description of our Gibbs-based algorithm developed . 
to partition hyperspectral imagery in an unsupervised manner. A 
description of the experiments using HYDICE imagery is given 
for both validating the algorithm and coniputing the spectral 
mixture process conditioned by the partition pirocess. 

II, Mathematical Framework 

A Bayesian model is introduced in this section to account for 

spectral mixing and the partitioning of homogeneous regions in 
hyperspectral imagery. The approach is model-based. Multiple 
conditional probability structures are used to bind together a 
hierarchy of underlying and unobserved processes that combine 
to generate a given scene. We assume that there is an underlying 
hierarchy of two processes 

X = {X';X''). (2) 

The process X^, called the partition process, will identify 
regions that shall be treated as homogeneous. The process X^, 
called the spectral mixing process (SMP), will model the re^ 
fleeted energy, possibly being that of a composite of materials. 

Then, what is actually observed is G, the observed hyperspec- 
tral image cube 

= + noise." (3) 

For each of G and X^\ there are two "coordinates," one spa- 
tial and one spectral. 

Notationally, we make the following definitions. Let Sj de- 
note the full-resolution image pixel lattice, i.e., the resolution at 
which the observations are obtained, which will also denote the 
full-resolution partition lattice Sp 

Sp = 5/ = < i < JVr^wb, 1 < J < ^pok} ' (4) 



where iVaows is the number of rows in an image, and iVcois is the 
number of columns. Let A denote the full-resolution wavelength 
lattice 

A = { A/ : Ai e (400 nm, 2500 nm) ,1^1] iVbands } (5) 

where A^bands is the iiumber of spectral bands in the hyperspecr 
tral image cube. Hence, for example, the full-resolution coordi- 
nates for the observed hyperspectral G are 5/ x A:^ 6 5/, A € 
A. Therefore, G{3, A), is the observed intensity at that "site" and 
wavelength. 

Below, we present the constructions of X^\ X^ ^ and G, 
Many of the explicit details are given in Section III. 

A. Construction of the X'* Process 

As stated in Section I, our primary objective is to develop 
a method for image segmentation that utilizes both spatial and 
spectral information. Hence, the construction of X^ is fundar 
mental. Our approach for the spatial-rspectral partitioning will 
be multigrid, and thus we define a collection of suMattices of 
5/> for the partition process 

0<z< — — . — ^,0<j< > (6) 

a ^ ) 

where cr is a parameter called the "grid resolution," which deter- 
mines the spatial resolution (coarseness) of the labeling process 
at some level in the hierarchy. 

Below, the term "site" will refer to a generic element of one 
of the aforementioned lattices. At full-resolution it would be a 
pixel, but on a cruder scale, it will be a square block of (full-res- 
olution) pixels. The value of a is defined to be the number of 
pixels per side of each block of pixels. A site might be com- 
posed of one material, but it could possibly be composed of nu- 
merous materials, and we need. to presume that we will not know 
this, or the identities of the materials, ahead of time. Moreover, 
the above sets, 5/ and 5yf\ will have associated graph struc- 
tures which describe, for a given element, its neighboring ele- 
ments. These specified structures allow us to sequentially for- 
mulate well-defined probability models for X^' (as a Markov 
random field [MRF]) as well as for X^ and G. 

Notationally, if 5 is one of 5/ or : S — {^i, . . . , }, 
where M is the number of sites on the grid, then an associated 
neighborhood system Q = \Qs^s € 5} is specified (details in 
Section III) which is a set of lattice graph structures of neighbors 
such that s ^Qs, and s € Qr if and only if r ^Q^. Within this 
setting, our goal is reduced to finding the value of X^', which 
maximizes Pr{X^|G) for the observed G. 

X^ is a discrete labeling process that associates a label with 
each site and is treated as an MRF. The realizations of Xj* are 
xf € r, where F = {1,2,..., iVLabeis}- We are using the 
general notation of upper-case "X" being a random variable 
and lower-case "a;" being a realizable value. A^Labeis is the total 
number of labels assigned to the partitions. The configuration 
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space for partitioning is then determined by iVLabeis and a (the 
label resolution), and denoted by 

= {x^ - (xf , S e S^r^) \xi' € {1, ... , iVLabels}} - 

Sp^ is the set of sites used in the partition process at a resolution 
(7. The development for is similar to that of [33] with the 
important addition of a spectral dimension. 

We treat X^' a Markov random field (MRF) with respect to a 
graph {6ir\g^<^)} so that 

Pr(^^ = x^') > 0 for all e (8) 
Pv{X^^x^\X^ = xi:,r^s) 

^Pr(jCf =a;M' = <,r€ef^-)) (9) 

Where {X^' =^ x^"} rneans {X^^ " ^f, , • - - ^ x'J^ }, 

Q^'^"^^ is the entire neighborhood structure for the partition 
process at resolution cr, and Qs^'^^ is the specific neighborhood 
at site s for the partition process at resolution a. 

We exploit the Gibbs equivalence theorem offered in [34] 
and with an appropriately defined graph {S^p\Q^^^^^}^ we can 
adopt a Gibbs representation 

ri{X^'^ x^,G^ g) = ie-*^(*'''^> (10) 

where Z ~ ^^i^ ^~{^/t)u{x^ ,g) j^^own as the partition fiinc- 
tion of statistical physics (not to be confiised with the "partition 
process") and U (x^\ g) is the total energy function representing 
the sum of all the interaction energy of each site with its neigh- 
borhood. We denote g to be the spatial and spectral realizations 
of G, which is in contrast to the spatial-only realizations of the 
simpler G discussed in [34]. Specifically, g^ are the spectral ob- 
servations at site s, and g are the collections of these observa- 
tions for all.the sites. Equation (I Oj is the global distribution for 
the entire image. 

An appropriate localized form of 6/ at a given site 5, 
(x^\g), is of practical value. The advantage of the Gibbs 
form is that we can compute the maximum a posteriori (M AP) 
estimate that maximizes the probability Pr(a;^^|^) by iteratively 
sampling from the local Gibbs distribution pertaining to each 
site € 5/f \ We use 

= i-e-*i'-(-''.f) (11) 

whereZ, ==X;a:f€re~^^^^^^^^^^'''^'''^^ andf/3 .(a:^\^) is the en- 
ergy interaction of site s 6 Sp^ with the neighborhood Qi 
Quite conveniently, the computations of Zg are not necessary 
for purposes of obtaining a MAP estimate. 

The result of the partition process is a set of partitions (de- 
composing the set S'p ^) 

S^{H,,(7er,5,c5jf>} (12) 



where Eq denotes the collection of all s in Sp, for which 
random variable Xf takes the value q. For each there will 
be a set, their union forming a decomposition of Sp . 

B. Maximization of L{/3\X^\X^) 
We define 

if-{/i(^>(A),...,^(^-^^^(A),A6A} - (13) 

as the set of material spectra, and C H as the set of material 
spectra at pixel site s. At a site s, h^^^ = {h^^\Xi)yXi e A) 
is the spectra for the A:th endmember if h^^^ € HsAn the case 
where there is one partition, Hg equals H, However, in the case 
where there are multiple jjiartitions, Hs may be strictly contained 
in if. JVends is the total number of fundamental materials. Recall 
that, unless we have accurate ground-truth information about the 
materials in a scene, the task of determining this number is not 
trivial, as discussed in Section I-B. The goal is to extract the ftin- 
damental materials, called endmembers, and associated coriipp^ 
sitions for materials imaged in a scene that cannot be observed 
directly, resulting in improved accuracy over existing isolated 
full-pixel and mixed-pixel methods. 

X^ is defined to be a wavelength-dependent random process 
governed by the spectral mixing relation 

<W= Ji''^-h^''H^) + nsW. (14) 

is the vector of proportions at a site 
s on a label lattice. The function rys(A) is associated with error 
due to variance in the endmember spectra. are the esti- 
mates of 0a^^ obtained through sonie type of constrained least 
squares or other means. Because the hyperspectral cube C gen- 
erates discrete samples of the spectrum at specific wavelength 
intervals, we can approximate the SMP process-using standard 
linear modelirig matrix notation 

X^=H,I3,-^V., ^€5/. (15) 

The basic SMA in (1) is the special case of there being only 
one partition, i.e., (14), as is, without the modification for the 
partitions. 

The physics of the spectral niixing phenomenon implies that 
the estimate of ffs should incorporate two constraints: a strict 
positivity constraining /Si^^ > 0 for k: /i^^^ € and a sum 
to unity constraint Ylk:hi'')eH ^^^^ ~ ^' T^^^se are properties 
exhibited by compositional data. 

Therefore, care should really be taken to treat the as com- 
positional data. -Statistical tools for analyzing compositional 
data have been developed [36]; use of these tools is envisioned 
in a more sophisticated process than the current one, as will be 
discussed at the end of this section. 

C: Construction of the Vi{X^\G) 

Now suppose we are working with high-quality data and have 
taken proper precautions sp that we can impose the assump- 
tion of uncorrupted data such that (3) can be written simply as 
G = X^ . Roughly speaking, without loss of generality, we are 
thinking of all the noise as being contained in the r/ term of (14). 
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We can speak of the full data as being (2), i.e., {X^, X^'). What 
is observed is not the full data, but rather only X^An principle, 
one could obtain the marginal probability density Fi{X^ \(i) by/ 
integrating out X^ (from the joint density Pv{X^\X^\P) ). 
One could then maximize the likelihood L{j3\X^ ). 

Because of the enormous inherent complications of such a 
complex integration (over all possible realizable partitions), we 
use the alternative two-step procedure, motivated by the expec- 
tation maximization (EM) algorithm [37] 

Step 1 ) Calculate the MAP estimate of denoted as X^, 
(analagous to the expectation step) 

= argiiiax Pr(X^'|Ji:0- 0^) 

Step 2) Calculate the maxinlum likelihogd estimate (MLE) 
of /?, denoted as p, (the maxirnization step) 

^^argmax (17) 

using the MAP estiniate of X^* as if it were "ob- 
served." ^ 
This is analogous to the EM algorithm, except since 

Pi:{X^\X^ ) does not depend on in pur case, we do not need 

to iterate as in the usual EM procedure. 

p. Extension of the Model 

Ultimately, a more sophisticated process is envisioned, 
namely 

X^iX^,X^,X^,X^) (18) 

where the components of X are a hierarchy of MRFs. X^ would 
be a stochastic line process identifying a boundary that is de- 
fined as an MRF with respect to a graph structure that defines 
its neighborhood system. Such line elements roughly "form" the 
boundaries separating the partitioned regions. X^ would be a 
stochastic compositional process that further constrains the end- 
member combinations and governs a relaxation of the mixing 
proportions /3s^^ within local neighborhoods. There is a rich 
statistical subject area concerning compositional analysis that 
can be used to develop techniques for determining this compo- 
sitional process [36]. The line process is envisioned to provide 
the capability to sever neighborhood structures between sites be- 
longing to different compositional phenomenon. Our goal is for- 
malized as the ability to simulate from Fr(Jf — x\G = g), via 
probabilistic relaxations in stages, i.e., {X^\G), . . etc., with 
the ultimate goal of replacing the process in (2) with the more 
sophisticated process in (18). The end result is a restoration of 
the hyperspectral cube froni which the local materials and their 
conipositions can be easily extracted. 

III. Description OF Algorithms 
A. Partition Algorithm 

The Gibbs-based algorithm implemented here is based on the 
framework just discussed. However, before one can implement 
the algorithm, a few practical matters have to be considered. 
Specifically in the following A I - AS, we define the energy func- 
tion of the Gibbs distribution introduced in (II), the associated 



neighborhood structures Qs, and the method of computation to 
. maximize the probability in (10). 
. The algorithm was designed to respond to dilferences in the 
shapes of spectral curves and to not be unduly sensitive to illu- 
mination differences. Consequently, the energy function below, 
measures the disparities in spectral angle between sites and their 
neighbors. This is in contrast to the traditional clustering algo- 
rithms, such as KMEANS and ISODATA, that use the Euclidean 
distance metric as a disparity metric (as well as a significantly 
different approach). 

I) Energy Function: The energy function used is 

{x'\g) = 6 (xf . (19) 



where Or^s = cos 




The delta function 6{xi,Xj) = 1 forxi ~ Xj. SLnd S{xi,Xj) = 0 
for Xi ^ . The quantity ^,-,3 is the spectral angle between 
s and r. If By, a < ^tiiresh, then SiQr^a) = —1, otherwise 
/(ft-.a) =^ Ij where the threshold A^thresh deterniines the sen- 
sitivity of the algorithm to differences in spectral angle. The ef- 
fect of this energy function is to accurnulate the value of /(^r,a) 
whenever the label of a site s is equal to its neighbor at r. 

The value of the threshold i^tiiresh is determined through ob- 
servation, specifically, it is determined by observing the differ- 
ences in spectral angle for samples of the sanie material versus 
the differences in spectral angle for samples of different mate- 
rials. If the imagery to be processed is calibrated to reflectance, 
then the value attained should be fairly universal, presuming 
a representative set of materials has been analyzed. If the im- 
agery has not been calibrated, then this value will depend on the 
range and scaling of the uncalibrated data. A proposed value for 
calibrated data is given in Section V based on observations of 
the spectral angles between the numerous spectra in a reference 
spectral database. 

Conceptually, one can think of the energy function as an 
accumulation of attractive and repulsive interactions between 
sites and neighbors. Whenever a site and its neighbor are in a 
common state (i.e., they are labeled the same), they experience 
an attraction for < /^thresh and they experience a repulsion 
for Br,9 > -ft^thresh. The attractive interaction contributes neg- 
ative energy and the repulsive interaction contributes positive 
energy. If a site and its neighbor are not in a comnion state, 
they do not interact, 

2) Neighborhood System: As indicated from the be- 
ginning, we implement a hierarchical multigrid approach. 
Therefore, the neighborhood systems change according to 
the grid resolution a. The algorithm proceeds in a sequence 
<^ = (^maxj ... ,8,4, 2, 1. However, throughout the sampling, 
we use a predetermined neighborhood structure that satisfies 
the condition that s e Oa i^nd s ^ Qr if and only if ?• 6 6s. 
The structure of the neighborhood system Qa forms a rather 
simple pattern. A neighborhood of a site will consist of near, 
intermediate, and far neighbors, and they are all at specific 
multiples of a. Our local neighborhood is either the four closest 
(above, below, left, and right), or eight closest (the surrounding 
perimeter) sites enclosing the site 5. The intermediate and far 
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neighbors form a similar pattern as the near neighbors but 
are at a further distance that is a specific multiple of cr. The 
intermediate and far neighbors are used to aid in obtaining a 
faster global solution and to enable the algorithm to remember 
label assignments computed by the algorithm at prior stages, 

3) ' Calculating the MAP Estimate: A simulated annealing 
technique is used here to obtain the MAP estimate by the 
Metropolis algorithm as described in [32], [34], and [38] with 
an eventual "fast freeze" at low temperature. For each iteration, 
we cycle through all the sites in an image. For a particular site 
s and iteration k, we compute 



7«(A;) =exp 



(21) 



The quantity 7a (A;) is the ratio of two Gibbs distributions, 
specifically, the ratio of the Gibbs distribution for a proposed 
configuration x^' and the Gibbs distribution of the current con- 
figuration where the cancels. The proposed configurar 
tion is identical to that of x^\ except for a proposal label 
change at site a*. If 7a (A;) > 1, then a transition is always made 
(the configuration of at iteration A: -f- 1 will be set to the 
proposed configuration). If 73(^') < 1, then the transition is 
made with a probability of 7s (/c). This event occurs whenever 
75 (A:) > ^, where ^ is a threshold probability that is chosen uni- 
formly in the interval (0,1), Otherwise, no transition occurs (the 
configuration of at iteration A: + 1 will remain unchanged). 
The temperature tk defines an annealing schedule and we use 
the one proposed in [32] and [34], namely 



ln(Hh/c)' 



C is a parameter corresponding to the maximum depth of a 
"well" in the energy function (19), which is discussed further 
in [32]. The value of C affects the rate of convergence. In our 
case, an upper boiind on this value can be computed as the value 
of (19), for which the spectral angle of the site-neighbor pairs 
are all different, but the labels are assigned the same value. A 
large value of C starts the computations at a high temperature, 
and it is theoretically appealing because a suiTiciently high tem- 
perature guarantees convergence to the true Gibbs distribution 
(a global solution). Computationally, a large C is undesirable 
because it results in very slow convergence and long computa- 
tion times. Selecting a small value of C tends to cool the system 
too quickly and causes the algorithm to get stuck in a local min- 
imum, where convergence to the true Gibbs distribution is not 
achieved. Because we are using a multigrid approach, the com- 
puting time can be reduced considerably by using larger values 
for C at the coarser grids, but smaller values at the finer grids. 
Practical values of C are discussed in Sections IV and V. 

B. Spectral Mixture Process Conditioned by Partitioning 
(SMPCP) 

Basic unconditioned spectral mixture analysis is the special 
case where the scene is treated as having one partition. In this 
one partition, there is one set of endmember spectra so that 



Hs =^ H, 3 ^ Sj. There are obvious disadvantages to this ap- 
proach. First, the estimates for (15), provided by a least 
squares form of (17), constrained or unconstrained, are most 
stable when there are only a small number of endmembers in 
the model. For typical scenes, there are often many. Using a 
small number of endmember spectra for such scenes results in 
large residual errors (poor estimates) for many sites, whereas 
using the larger ftill set of endmember spectra results in unstable 
estimates. Another disadvantage is that there are certain bland 
spectra, such as asphalt, water, shade, and metal that are very 
difficult to unravel in a mixture model, and lumping them to- 
gether in a model results in very poor results. 

With multiple partitions we can assign si different set of 
endmember spectra to each partition within the framework 
described in Section II-B. This method will be called "spectral 
mixture process conditioned by partitioning" (SMPCP), and it 
is consistent with the framework defined by (12)-(17). Ideally, 
we would like to determine these different endmember spectra 
automatically through methods such as [24]-[29]. However, as 
mentioned in the introduction, this is a very ambitious task. In 
an effort to control the sources of error, we extract image^de- 
rived endmember spectra, and based on kiiown ground-truth of 
the scene, manually assign these endmembers to partitions. 

To accomplish the conditioning, partition masks are 
constructed using the output class maps generated by the 
Gibbs-based algorithm.' A set of permissible endmember 
spectra is assigned to each partition to obtain Estimates 
of pi^^ in (15) are computed using constrained spectral 
mixing, where the masks are used to determine which set of 
endmember spectra are permissible in the model. The spectral 
toolbox available in ENVI [39], which applies the sum to unity 
(22) constraint, is used to compute the Pa^^ estimates on the masked 
data in a sequence of mutually exclusive operations. 

IV. Description of Experiment 

A. Data 

Two H YD IC E scenes extracted from a data run collected 
at 10000 feet altitude near Fort Hood, TX, on May 16, 
1996, are used in the experiment. Each scene consists of 
300 samples x 320 lines of pixel vector data with a spatial 
resolution of approximately 2 m. HYDICE Band 49 of Scenes 
1 and 2 are shown in Figs. 1 and 2, respectively. Complete 
HYDICE hyperspectral image cubes have 210 bands in the 
approximate 400-2500 nm spectral range. However, because of 
the overwhelming degradations that occur in the atmospheric 
windows and occasion line degradations occurring in certain 
bands, only 1 1 7 bands were actually used in the analysis. Fig. 3 
shows the regions (shaded) that were excluded from the pro- 
cessing, as well as grass and soil spectra extracted from Scene 
2. The 1 17 bands of data actually used in the computations, for 
each scene, were calibrated to reflectance by the empirical line 
method based on the known spectra of the calibration panels 
that are evident in Scene 2. 

Ground truth, documenting the identity of some of the ground 
features in Scene 1, is based on data collected during a site visit 
to the area during April 2000. As part of our experiment, the 
data from this site visit (a separate effort) were compiled as a 
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Fig. I . Scene 1, HYDICE Band 49. 




Fig. 2. Scene 2. HYDICE Band 49. 




bf^<«if ad [ti l I 



Fig. 3. Spectral signatures of grass (Ul) and soil (#2). 

graphic overlay shown in Fig. 4, with an associated descriptive 
table, Table I. Because of the time difference between scene ac- 
quisition and the site visit, not all portions of the scene could 
be documented with certainty. In addition to the ground-fea- 
ture data, spectral data of a diverse range of materials, hiea- 
sured by field spectrometers and collected independently of this 




Fig. 4. Ground truth overlay of Scene 1. 



. TABLE I 

Description of the Ground Truth (GT), as Displayed in Fig. 4 

gTMW DeicriottoncttheAirm 

A AsphaJt parking lol, newly paved! 

B Rooftop of a dq)aitnKnt store Jight>gjsy asphalt shingjc^^ 

C Concrete pavement to loading dock; some tire inaika. ' 

D AsjAalt parte ing lot, 

E , Church and school uith metal rooftop: 

F Private residettcc, slaiid-alone. 

G ' School. ropfU>p coated with bubbly light-gray maierial, top edge is grawl. 

H Apvtmmt fauldings with asphalt^ 

i Cornig^ed steel buldings. 

J Shade of building. 

K Asphalt parking lot surrounding an apartment coniplex. ' 

L^. Church with dark brown asphalt-shingled rooftop. ' . 

M Tenniscourt, concrete with blue/^oy coating, and lots of dirt. 

N ^ Playground area with exposed soil. ' , 

O Encircted area ifi reudentiet with inoderately-priced homes and concrete driveways. 

P A^alt cpad netwock within cndrckd area. 

Q LighMiiiear pattern is eiqNMed bedrock. 

R Tnees - Deciduous holly, Bois D Arc. Osage Orange. Texas Oak. and Bow Wood 

S- Ti>ccs. Red Bud, Texas Oak, and Red Bud 

t . Scattered junipers throughout the grassy areas. 

U Tall wild grass. . ' - 

V A^alt road intersection with exposed soil along the shoulders 
W Playground with exposed soil, gravel, and wmc asphalt. 

X Tall wild gnus. * 

Y Exppttd soil ffid bediodc. 



experiment at other geographic locations, were used to deter- 
mine an appropriate quantization threshold for the spectral angle 
disparity metric, as discussed in Section IV-B. Unfortunately, . 
groundrmeasured spectral data within the study site were tiot 
available. 

Of the two scenes, Scene 1 has the type of complexity where 
it is more appropriate to apply the entire Gibbs-based condi- 
tioned spectral mixture analysis^ Due to the nature of Scene 
2 (almost exclusively vegetation), this scene does not stand to 
benefit much from a spectral mixing process conditioned by 
multiple partitions. In the context of conditioning the spectral 
mixing process, one partition is sufHcient. Scene 2 is really an 
exaniiple of the kind of scene that is appropriate for either tra- 
ditional SMA or traditional full-pixel methods. In our experi- 
ments. Scene I was used to test the entire process. Scene 2 was 
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TABLE n , 

List of Endmember Classes Used in the Spectral Mixing Trials 



Fig. 5. ISODATAcluster map of Scene l . ' 

used primarily to test the partition algorithm in a different scene 
environment. 

B. Characterization and Validation of the Partitioning 
Algorithm 

Computer code was written in C++ to implement the Gibbs- 
based algorithm in the manner described in Sections II and III. 
The spectral angle disparity metric was implemented as defined 
in (19) and (20). In order to determine an appropriate value of 
the quantization threshold A'tUresh, the spectral angles between 
ail the pairwise combinations of the spectra in the reference 
spectral library, just mentioned previously, were coniputed. 

The two scenes were processed in multigrid sequences of 
runs. Because the algorithm is new, some experimentation with 
the parameters is necessary to get a flavor of the algorithm's . 
behavior with different settings. The parameters that affect this 
behavior are a (grid resolution), C (controls the annealing 
schedule through the temperature tk\ i^totai (final number 
of iterations through the image), and iViabeis (the number of 
• partition labels). After sonie initial probing of the parameters, 
three trials .(Trials A,. B, and C) were perfomied on Scene 1 
and two trials (Trials D and E) were performed on Scene 2. 
The results are documented in Section V, >yhich include the 
computation times for the runs. 

The unsupervised ISODATA clustering method is used as a 
reference method, and serves to provide motivatiori for using the 
Gibbs-based approach. Fig. 5 shows the results of this algorithm 
on Scene 1, for 6 classes. • 

C. Spectral Mixing Experiment 

Two spectral mixing trials were performed using Scene 1, 
which is considered the more appropriate of the two scenes for 
the conditioned spectral mixture analysis. In the first trial, con- 
ventional spectral mixture analysis was performed in the con- 
ventional manner with the ENVl spectral toolbox. Ground truth 
was used to determine appropriate endmembers. The spectral 
data of the endniembers were extracted from the scene, based 
on known locations of the endmembers. The class names are 
listed in the first column of Table II. 



NoPwtlttonlng 

Bright Roof (B) 
Concrete/ Exposed Soil (Y) 
Asphalt Pavement (A) 
Asphalt Shingles (H) 
Trees (R) 
. Grass (X) 
Shade (J) 



Bright Roof (B) 
Concietc/ Exposed Soil (Y) 
Asphalt Pavement (A) 
Asphalt Shingles (H) 
Shade (J) 
Qt?ss(X) 



Partition B 

Concrete/ Exposed Soil ( Y) 
Trees (R) 
Grass (X) 
Shade (J) 



■ „ TABLE III 
Control Parameters for the Partitioning Trials. The Control 
Parameters for the Scene I Trials (Trials A, B, and C) and the Scene 

2 Trials (Trials D and E) are Shown Below. For Each Grid 
Resolution a. We List the Parameter C, the Number of Annealing 
Iterations, the Neighborhood Structure (the Pixel Distance of the 
Neighbors from a Site), and the Running Time on a 500 MHz Laptop 
Computer, (a) List of Control Parameters for Trial a. (b) list of 
Control Parameters for Trial B. (c) List of Control Parameters for 
Trials C and p, (p) List of Control Parameters for Trial E 
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In the second mixing trial, a simplified approach of condi- 
tioning the spectral mixing process on the partition process was 
applied. Ground truth was used to aid in the assignment of end- 
members to the partitions that were generated in Trial A by the 
Gibbs-based algorithm. Based on scene knowledge and the par- 
tition results, the six partitions resulting from Trial A were col- 
lapsed into two partitions (A and B), representing two types 
of phenomenon: natural/vegetation and culturaUbuilding ma- 
terials. The endmember spectra extracted in the first mixing 
trial were reassigned so as to belong to Partition A or Parti- : 
tipn B. These assignments are listed in the right two columns 
of Table II. Note that although the partitions are mutually ex- 
clusive, the class assignments often overlap. Also, note that the 
letters in parentheses of Table II indicate the regions of the scene 
where the spectra were extracted, as shown in Fig. 4 and Table I. 

V. Results 

Fig. 5 and Table IV(a) show the results of the ISODATA clus- 
tering method applied to Scene 1 for six classes, using the same 
bands selected for the Gibbs-based sequences discussed below. 
Although ISODATA provides a fairly good partitioning of some 
materials in the scene, it has a very bad problem confiising 
asphalt with significant portions of vegetation in the forested 
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(b) 



(e) 




(d) 



(0 



Fig. 6. Partitioning results for spectral angle disparity, Scene 1 . (a) Random initialization a = 8, (b) results for trial A cr = 8, (c) results for trial A a 
results for trial A <r == 1 . (e) results for trial A <y = 1, (£) results for trial B a = 8, and (0 results for trial Ca = 1, 



4.(d) 



area, denoted by "R" in the ground-truth map (Fig. 4), and has 
difficulty distinguishing between grass and the trees in "R." 
These problems are seen explicitly in Table IV(a). Out of the 
3019 pixels extracted frpm the "R" region, 580 were placed 
into ISO-Ol, 1426 were placed into ISG-02, and 1012 pixels 
were placed into ISO-04. The cluster associated with ISO-04 
can be associated with asphalt by observing Table IV(a), Row 
A. The cluster label ISO^02 can be associated with Grass by 
observing Table IV(a), Row U. For purposes of the subsequent 
conditioned spectral mixture process, confusion between grass 



and trees is not problematic. However, confusion between as- 
phalt and (slightly shaded) vegetation is very problematic. The 
clustering also does not provide a smooth mapping of regions, 
spatially, and exhibits a significant amount of speckling. 
' Based on observations of spectral angle distances between 
the various pairs of materials, we set /^thresh — HA which re- 
mained constant throughout the partition trials. In most of the 
pairwise signature combinations, there was an obvious gap be- 
tween similar materials and different materials, except for com- 
binations involving a vegetation signature. Unfortunately, veg- 
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■■' . (a) 

Fig. 7. Partitioning results for spectral angle disparity. Scene 2. (a). results for 

elation has a fairly large variance which rnakes the threshold 
boundary a bit fuzzy. The exact value was observed to not be 
terribly critical. However, valiies smaller than iCtiiresu ~, 8.0 
would make the algorithni overly sensitive to feature diflFerr 
ences, and values larger than A^thresh =^ 15.0 would make the 
algbrithm insensitive to feature differences. 

The partitioning algorithm was tested for A^iabeb = 6 and 
-/Viabeia = 8 partition labcls, with no significant difference in. 
partitioning capability. Although ^labels = 6 is obviously too 
small to account for all the possible materials in a truly global 
labeling of the scene, it is also obvious that we are unlikely to 
determine the true number exactly, in an unsupervised mode. 
This number is usually sufficient, however, to account for any. 
local variations of materials in a geographic area, and, in the 
case of moderately complex scenes, to capture the predominant 
materials globally. We should emphasize that the goal of the par- 
titioning algorithm is not to determine the identity of partitions: 
The identification of the composition in a scene occurs through 
the conditioned spectral mixing process. Therefore, determining 
the exact number of partitions is not critical. However, if the 
algorithm is applied in a supervised rrianner (an alternative ap- 
proach that was suggested in Section I), the number of partitions 
could be increased or decreased according to scene knowledge. 

Table III lists the values of parameters used to generate the 
partition maps shown in Figs. 6 arid 7, as well as the computation 
time that is required on a 500 MHz G3 -based laptop computer. 
The range C (1.0^16.0) was investigated, which is below 
the upper bound discussed in Section III-A.3. The specific value 
affects the rate of convergence. As an example, for C = 16, the 
destiny of the partitioning results seemed to be controlled by the 
first 1 50 000 iterations, and for a value for C = 2, the destiny of 
the partitioning results seems set within the first 2000 iterations. 

Fig. 6 shows the results of the Gibbs algorithm applied to 
Scene I. As indicated in Table III, each trial is executed as a 
multigrid sequence of runs, proceeding from a course resolutipn 
to finer and finer resolutions. The output of one stage is used 
tp initialize the next stage at a finer resolution. Fig. 6(a) shows 
the random initialization of Trial A. Figs. 6(b>-{d) are the re- 
sult of Trial A using the parameter values listed in Table Ill(a). 
Summing the corresponding run times, the total computing time 




(b) / 

trial D a = 1 , and (b) results for trial E<r±:l. 

TABLE . IV . • ' 

Confusion Matrices of Results. Labels in the 15 Rows, are the GT 
Labels, and Labels in the Six Columns are the Partitioned Labels, 
As an Example of How to Read Table IV (a), the Region of Groljnd , 
Truth Extracted for L Consisted of 82 Pixels. Of the 82 Pixels, 44 

Pixels were Assigned to Cluster Label lsp-03, 37 Pixels were 
Assigned to Cluster Label lSO-04, and One Pixel was Assigned to 
the Cluster Label ISO-06. (a) Results for the ISODATA Algorithm 
and (b) Results for the Gibbs-Based algorithm: Trial A 
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of the multigrid sequence for Trial A was quite long (8 h, 45 
min). Fig. 6(e) is the final result of Trial B using the param- 
eter values in Table 111(b), A comparison of Fig. 6(e) and (d) 
shows that the result of Trial B is almost identical to the result 
of Trial A. The total sequence was computed in far less time 
(58 min), with the difference in execution between the two runs 
occurring only at the last two stages. Fig. 6(f) is the final re- 
sult of Trial C, achieved in even less time (28 min) using the 
values in Table III(c). A comparison between Figs. 6(f) and (e) 
shows some labeling differences. However, the structural infor- 
mation provided by the partitions that will subsequently be used 
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Fig. 8. Abundance maps for bright rooftops, (a) Basic SMA (X'^jG). and (b) SMPCP (X^IX^'G). 



.(b) 




(a) 

Fig. 9. Abundance maps for trees, (a) Basic SMA (A''' |G); and (b) SMPCP {X^X^'G). 



(b) 



by the conditioned spectral mixing process is not significantly 
different. 

Table IV{b) shows quantitative results for Trial A. A com- 
parison of Tables iV(a) and (b) provides justification that the 
Gibbs-based partitioning gives advantages with respect to the 
partition map extracted with ISODATA for the subsequent con- 
ditioned spectral mixture process. Observing rows A, C, and 
P [Table IV(b)], the label LBL-04 can be associated with as- 
phalt or concrete. Observing row U, LBL-01 can be associated 
with grass. Out of the 3019 pixels extracted from the "R" region, 
2800 pixels are assigned a single label (LBL-02), only 17 pixels 
are labeled as LBL-04, and 152 pixels are labeled as LBL-01. 
The troublesome confusion between (slightly shaded) vegeta- 
tion and asphalt that plagues ISODATA is not a problem for the 
Gibbs algorithm. 

Nevertheless, there are still some problems in distinguishing 
certain urban categories. For example, LBL-03 is assigned to 
both "H" asphalt-shingle rooftop and "K" asphalt parking lot. 



Also, a problem in the global labeling of small and/or skinny 
features that can be observed in Fig. 6 is not represented in this 
table. These problems will be addressed in future work. 

Fig. 7 shows the results of the Gibbs algorithm for two trials 
(Trial D and Trial E) applied to Scene 2. Fig. 7(a) is the result of 
Trial P, running the algorithm with the same parameter values 
as those used in Trial C. Fig. 7(b) is the result of Trial E, run- 
ning the algorithm with the values in Table Ill(d). A compar- 
ison of Figs. 7(a) and (b) shows no significant difference in the 
structural informatipn provided by the partitions and only a few 
labeling differences. The total run time for Trial E was 10 min 
and 39 s. 

From a theoretical viewpoint, A^totai should be very large to 
guarantee a global solution. However, from a practical view- 
point, we would like to reduce the value as much as possible, 
particularly when the algorithm is operating at a fine resolution 
{a = 2, 1). Also, from a practical viewpoint, we should keep the 
starting teniperature, controlled by the parameter C, somewhat 
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(a) . = 

Fig. 10. Abundance maps for grass, (a) Basic SMA (X^jG), and (b) SMPCP (X''|X''C?), 



(b) 





(a) . 

Fig. 1 i. Abundance maps for asphalt, (a) Basic SMA {X'^lG), and (b) SMPCP (X' |X^G), 



cool at the finer resolutions for two reasons: I) to reduce the 
number of iterations and 2) so that the algprithm retains much of 
the information provided by the coarser resolutions. According 
to the trials discussed above. Tables III(c) and (d) show a list of 
parameters that provide practical computing tinies without ad- 
versely affecting the final solution. 

Overall, the Gibbs-based algorithm generates partitions that 
accurately represent structured regions in the image. For all of 
the bigger regions, pixels with the sanie label represent the same 
type of phenomenon on the ground, globally across the image, 
and regions with different labels correspond to different phe- 
nomenon. However, sonie problems can be seen in a few of 
the smaller regions. Although the contiguous pixels in smaller 
partitions represent the s£une phenomenon, locally, some small 
partitions, separated by a certain amount of distance and given 
the same names, actually represent different phenoniehon. Con- 
versely, spfiie sniall partitions, separated by a certain amount of 
distance and given different names, actually represent similar 
phenomenon. According to ground truth (Fig. 4), the algorithm 



was not able to distinguish the apartment rooftops (H) from 
the surrounding asphalt parking lot (K). However, note that the 
rooftop material also contained asphalt (asphalt shitigles). We 
could set the algorithm to have greater sensitivity by loweriiig 
the value of the threshold i^thxesh- However, this would have 
the undesirable consequence of partitioning the scene into too 
many regions (particularly, the vegetation), for purposes of the 
conditioned spectral mixture process. 

The problem of global labeling is not unexpected because 
the Gibbs algorithm, by construction, tends to emphasize local 
consistency, spatially, by virtue of the MRF formalism. This 
problem is not as prevalent in algorithms such as KMEANS 
and ISODATA. However, these other algorithms suffer in an- 
other important way because local consistency is ignored. From, 
a perspective of the proposed Bayesian model, the behavior of 
the Gibbs-based algorithm is preferred, because the ultimate 
identification of materials in the scene is determined by other 
model components. Also, recall that we are investigating this 
approach using a random initialization scheme (in some sense. 
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> a "worst-case scenario," as mentioned in Section I-C). Alters 
native initialization schemes should substantially improve the 
global labeling. 

The fractional abundance maps on the left side in Figs. 8-1 1 
show the results of performing SMA using one partition with the 
endmember classes listed in Table Il^Column A. The fractional 
abundance maps on the right side in Figs. 8-1 1 show the results 
of performing SMPCP, as described in Section IV-C, using the 

. endmeniber classes listed in Table II, Columns B and C. Bright 
areas correspond to an estimated high fraction of the material 
being represented in the abundance image. Figs. 8^1 1 are the 
abundance maps for Bright Rooftops, Trees, Grass, and Asphalt, 
respectively. 

The results show an obvious iniprovement when using the 
SMPCP method. The most noticeable improvement is how the 
unwanted asphalt abundance, prevalent in Fig. 1 1 (a), is reduced 
by the SMPCP in Fig. 1 1(b). In the past, the poor estimation of 
asphalt abundance by basic SMA in regions of shade, vegeta- 
tion, and water has been a notorious problem. SMPCP ofiers a 
practical way to solve this problem. Overall, the SMPCP abun- 
dance maps in Figs. 8-1 1 are sharper, cleaner, and place better 
emphasis on the material of interest for any specific map. 

VI. CONCLUSIONS 

The Gibbs-algorithm, using spectral angle as a disparity 
metric, partitions the imagery into nicely structured regions^ 
that are appropriate for input into other subsequent prpcesses in 
the proposed Bayesian model. The partitions represent homor 
geneous regions without being overly sensitive to illumination 
changes on objects in the scene. There were some problems 
in achieving a global labeling for some of the smaller regions, 
which should be solvable by using alternative initialization 
schemes, rather than using random initialization. This problem 
did not adversely affect the ultimate performance of the 
combined process. 

An improved spectral mixture process was demonstrated 
using the SMPCP. Partitioning the data provides a mechanism 
for reducing the dimensionality of the SMP at each site in 
an image that provides for more stable abundance estimates. 
Restricting the allowable endmembers in a partition also elimi- 
nates false predictions, particularly for endmetnb'er spectra that 
are bland in character, such as asphalt. 

The success in applying the Gibb-based and SMPCP algo- 
rithms validates the concept of using a Bayesian model and a 
hierarchy of structure-inducing processes and provides the mo^ 
tivation for developing more sophisticated Bayesian approaches 
as proposed in Section II-D. 
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Abstract - An unsupervised Gihbs-based simulated 
annealing approach to partitioning hyperspectral imagery 
into homogeneous regions is investigated, where spatial 
consistency is imposed on the spectral content of sites in 
each partition. Maximum A Posteriori (MAP) estimates 
are computed through the use of a Gibbs distribution 
defined over a neighborhood system. The algorithm is 
investigated as a multigrid implementation that uses a 
non-traditional extended neighborhood system to improve 
global labeling and reduce computational intensity. An 
energy function is introduced which models spectral 
disparities in an image using a combined spectral angle 
and Euclidean distance measure thai is defined with 
respect to the neighborhood system, the algorithm is 
focused on terrain mapping applications using 
hyperspectral imagery containing narrow bands 
throughout the 400-2500 nanometer spectral region 
(visible through short-wave infrared). 

Keywords: Bayesian, Markov Random Field, partitioning, 
clustering, hyperspectral, annealing. 

1. INTROPUCTION 

We investigate an unsupervised method of segmenting 
hyperspectral imagery into regions that we shall treat as 
homogeneous. Qur goal is to smoothly label picture, 
elements (pixels) in a scene according to their material 
type without the use of a reference spectral database or the 
need for supervised training. Our approach will use the 
paradigm of Bayesian image analysis. 

A mathematical development of the underlying theory, 
applying the Bayesian framework to Random Processes, is 
presented in [1]. Bayesian approaches using Markov 
Random Field (MRF) models have provided algorithms for 
producing unsupervised partition maps and detecting 
boundaries on imagery containing one band [2]. Geman et 
al [3] use the Gibbs sampler and Metropolis methods in 
their algorithms for the restoration of imagery, 
partititioning, and boundary detection, to compute 
maximum a posteriori (MAP) estimates. A method to 
increase the efficiency of global optimization within a 
general MRF framework for image processing problems 
has been developed, the Renormalization Group algorithm 
(RGA), that performs a hierarchical, multiscale, coarse-to- 
fme analysis of imagery, for the restoration of degraded 
images [4]. Although multiband imagery has been treated 
to spnie extent iti the Bayesian comniunity, the 



significantly different (spectral) character of hyperspectral 
imagery has not really been addressed, 

In the present effort, we discuss a multigrid Gibbs-based 
simulated annealing algorithm for. partitioning 
hyperspectral irnagery. This effort can be considered an 
extension, spectrally, of the spatial approach used in [2,3], 
with additional increased efficiency provided by a 
multigrid implementation that uses a special case of the 
general course-to-fine grid structure used in the RGA 
approach [4]. A non-traditional extended MRF 
neighborhood system is introduced for improved global 
labeling. 



II. Description of Algorithm 

A, Framework 

A Bayesian model provides a framework to construct the 
partitioning of homogeneous regions in hyperspectral 
imagery. The process , called the "partition process," 
will identify regions that shall be treated as 
"homogeneous," in the sense that each region contains 
some culturally similar phenomenon (not necessarily a 
single material type). What is actually observed is G, the 
observed hyperspectral image cube. For G there are both 
spatial and spectral coordinates. 

Our approach for the partitioning is multigrid, and thus we 
define a collection of lattices that provides a finite 
sequence of coarser and coarser grids: 

S^^O _^ ^ _^ s^^'-^\ where . 



4^'^={(/^. + Uy^r, + l): 



o</$%^,o< y<^%^} (1) 



N^^q is the number of stages in the sequence, CT„ is a 
parameter called the "grid resolution" which determines 
the spatial sampling of the algorithm and the spatial 
resolution (coarseness) of the labeling process at stage n 
(n = \,...,Ng^^) of the multigrid processing, Nf^^^^ is the 
number of rows in an image, and iVc^/^ is the number of 

columns. 5^^"^ is the set of sites that will be used in the 
partition process at a resolution <7„ . Our grid structure is a 
special case of the general course-to-fine grid structure 
used in the RGA approach to iniage processing [4], where 



the coarsest level contains all the large-scale features and 
each finer level contains correspondingly smaller-scale 
features. The term "site" refers to a generic element of a 
lattice, as specified in (1). At full-resolution a site 
corresponds to a pixel, but at a coarser scale it corresponds 
to a square block of pixels. The value of (T„ is defined to 
be the number of pixels per side of each block of pixels. 
At the finest resolution cXj = 1 . We denote the full- 
resolution image pixel lattice as Sj, i.e., the resolution at 

which the observations are qbtaiiied, where 5/ = 5?' 

Now if S is one of Sf"^: S ^ }, where M is the 

nuniber of sites on the grid, then an associated 
neighborhood system Q^^^ ^ - {§s,s e S} is specified 
(details in Section IIB) which is a set of lattice graph 
structures of neighbors, such that s Q^, and sgQ^ if and 
only if r e Q^. 

X^is a discrete labeling process which associates a label 
Nvith each site, and is treated as an MRF. The realizations 
of are jcf eP, where T ={l,2,...,iVi^f,^;Jis the set 
of partition labels. Ni^^^^^ is the total number of labels 
assigned to the partitions. . The configuration space for 
partitioning is then determined by Nj^fy^i^ and <t„ (the 
label resolution), and denoted by 

Q<;»>. {;c''=(xf, 5eS^')lx;e{l,...,iV^,,j}. The 
development for X^ is similar to that of [2] with the 
important addition of a spectral dimension. We treat X^ 
as a Markov Random Field (MRF) with respect to a graph 

. where Q^^^^ is the entire neighborhood 

structure for the partition process at resolution a„ , and 
is the specific neighborhood at site s for the partition 
process at resolution a„. We exploit the Gibbs 
equivalence theorem [3] which recognizes the MRF as a 

- random process with a Gibbs distribution. The advantage 
of adopting this approach is that we can now conipute the 
maximum a posteriori (MAP) estimate that maximizes 
the probability Pr(x'' I g) by iteratively sampling from the 
local Gibbs distribution pertaining to each site 5€5p^"\ 
We use 

^Pr(X;^xf IXf =;c;,reg,) 

(2) 

where {x'' =x''} means {x,^ = x,^,... -xf^}, 
Z, ^ X e^^^'^'^'^^K and U,(x^:g) is the energy 



interaction of site s^Sp" with the neighborhood Q^, 
Quite conveniently, the computations of are not 
necessary for purposes of obtaining a MAP estimate. We 
denote g to be the spatial and spectral realizations of G, 
which is in contrast to the spatial-only realizations of the 
simpler G discussed in [3]. Specifically, are the 
spectral observations at site s, and g are the collections of 
these observations for all the sites. 

.A simulated annealing technique is used to obtain the 
MAP estimate by tHe Metropolis algorithm as described in 
[1], [3] and [5]. In computing the MAP estimate, we 
investigate an approach that is consistent with aspects of 
the Renormalization. Group Algorithm (RG A) method 
developed for the restoration of degraded images [4]. The 
approach, called "unconstrained simulated annealing," 
uses a grid-coarsening and cascading technique that is 
consistent with the RGA approach. Specifically, the fine- 
to-coarse sequence, described above, with partition lattices 
defined by (1) is reversed to a coarse-to-fine sequence: 

S^p'''^'^ ^ ^ . . , 5{f 2) _^ ^y). The results of stage are 
used to initialize the algorithni at the next finer stage 
/." In the interest of achieving better globalized labeling, 
we investigate making more than one pass through such a 
rnultigrid sequence. A constrained implementation of the 
algorithm is also currently under develppmeiit. 

B, Neighborhood System 

The typical MRF neighborhood is compact (i.e., the closest 
four or eight neighbors to a site). However, in order to 
extend the reach of the neighborhood to encourage a nriore 
global labeling without a tremendous increase in 
computation, our algorithm uses a predetermined local, 
intermediate, and far neighborhood structure that still 
satisfies the (if and only if) neighborhood condition, 
mentioned above. Local, intennediate, and far are relative 
to a given scale (CF^). 

The structure of pur neighborhood system forms a 
simple symmetric pattern. The neighborhood of a site will 
consist of near, intermediate, and far neighbors, and they 
are all at specific multiples of a„ . Our local neighborhood 
is the eight closest (the surrounding perimeter) sites 
enclosing the site s. The intermediate and far neighbors 
form a similar pattern as the near neighbors but are at a 
further distance that is a specific multiple of cr„. The 
intermediate and far neighbors are used to aid in obtaining 
a faster global solution and to enable the algorithm to 
remember label assignrnents computed by the algorithm at 
prior stages. Of course, different neighborhood structures 
will produce, when the conditional probability measures 
(expression (2)) are pieced together, (potentially) different 



joint probability measures on the overall configuratipn 
space. 

C Energy Function 

The energy function, defined below, depends on difference 
measures, called "disparity measures," which indicate the 
dissimilarity between a site s and a neighbor r. Our energy 
function shall combine two disparity measures that 
measure spectral difference in two very different ways: the 
spectral angle and the Euclidean distance; 

For a pair of sites at the site locations r and s, consider a 
set p^^={Z>^J],Z>J^^^^}, where the elements of D^^ are the 
spectral angle Z>J]^ and the Euclidean distance D^^ 
between the' sites s and r. The individual distance 
measures pJ'J, 1=7,2, are defined as follows: 

P^=sgr{ ^'f-'i^-''' ) : ; (4) 

For the set D,^={D^},D^^}, we define the energy 
function 

l5(^;,^r) niax{/,(Djj>),D<S (5) 

where the delta function S(Xy,x^)=\ forA;^=jCj , and 
^(jc^,jc^)=0 for Xr'^Xg. The function //(0[']) maps 
to the set {-I, 1}, where /,(D^;])^-l indicates a 
non-significant distance between r and 5, and a value 
/^(dJ'^) = 1 indicates a significant distance. If 
D^;]<^^^ then /^(D^;>) = -l, and otherwise 
/.(Dj;]) = l, where the threshold K^^resh determines the 
sensitivity of the algorithm to the Z)^j] distance measure. 
The effect of this energy function is to accumulate the 
value of /, (D^']) whenever the label of a site s is equal to 
its neighbor at r. 

The value of the threshold Kj^hresh determined through 
observation; specifically, K^hresh *s determined by 
observing the differences in spectral angle for samples of 
the same material verses the differences in spectral angle 
for samples of different materials, and K^hresh is 
determined by observing the differences in Euclidean 
distance for samples of the same material verses the 
differences in Euclidean distance for samples of different 
materials. If the imagery is calibrated to reflectance, then 



the value attained should be a fairly universal value, 
presuming a representative set of materials has been 
analyzed. If the imagery has not been calibra.ted, then this 
value will depend on the range and scaling of the 
uncalibrated data. . 

This energy function is designed to respond to disparities 
as ' measured by the spectral angle measure and the 
Euclidean distance, Our reason for using a combined 
measure is that we seek to increase the sensitivity of the 
algorithni without introducing undesirable clutter. Rather 
than using a single measure with a very sensitive threshold 
^Thresh* ^^^^ would Tcspond to onc specific type of 

spectral difference, we use two nieasures that respond to 
different kinds of spectral difference with less sensitivity. . 

III. Validation 

A HYDICE scene is used for validation that was extracted 
from data collected over a geographic area near Fort Hood, 
Texas at a 10,000 feet altitude, on 16 May 1996. Figure 1 
shows the scene, consisting of 300 samples by 600 lines of 
pixel vector data with a spatial resolution of approximately 
1.5 meters. Because of the severe degradations that occur 
in atmospheric windows of the spectrum and sensor 
degradations occurring in certain bands, only 117 out of a 
possible 210 band were used to compute the disparities 
within the energy functions. These 1 17 bands of data were 
calibrated to reflectance by the Empirical Line Method 
based on the known spectra of the calibration panels that 
were located in a region outside this particular scene. 

Ground truth, describing certain ground features in the 
scene, was based on data collected during a site visit to the 
area during April 2000. In addition to the ground-feature 
data, spectral data of a diverse range of materials, 
measured by field spectrometers and collected 
independently of this experiment at other geographic 
locations, were used to determine an appropriate 
quantization threshold for the spectral angle and Euclidean 
distance disparity metrics, as discussed below. Ground- 
measured spectral data within the study site were not 
available. Figures 2 and 3 show the results of two trials: an 
ISODATA clustering run (control method) and a Gibbs 
unconstrained multigrid run, respectively. Both trials are 
set to generate six classes. The partitioning algorithm is 
initialized randomly, and it then proceeds unconstrained in 
a two-pass multigrid sequence with grid resolutions 
<7„ =16,8,4,2 16,8,4,2,1. That is, the algorithm 
proceeds down from the coarsest grid to a grid resolution 
of ^2 " 2, at which point it subsequently steps back again 
to the coarsest grid and then proceeds to (7^ = 1 . The 
distances of the intermediate and far neighbors correspond 
to fixed multiples of the grid resoliition, 4cr„ and 8<7„ , 



respectively. Table I provides commission errors as a very 
brief summary of the quantitative results: Due to space 
limitations the original confusion matrices are not shown, 
and the numerous ground truth (GT) classes are 
consolidated to three classes. . 




Figure 2. ISODATA Clustering - 6 classes. 




Fi^re 3« Gibbs partitioning - 6 classes. 



Table L Errors of commission. 



Commission Errors 



Trial 


Trees 


Grass 


Asphalt 


ISODATA 


22.48% 


1.32%" 


68.66% 


Gibbs 


2.60% . 


2.79% 


24,61% 



Figure 3 shows a partitioning of the scene that has a 
significantly smoother appearance, with more well-defined 
structures, as compared to Figure 2. Quantitatively, Table 
I shows the Gibbs algorithm provides a significant 
decrease in the commission errors for Trees and Asphalt. 



One negative point is that the Gibbs algorithm (initialized 
randomly) appears to have a problem labeling small skinny 
features. , 

IV. Conclusions 

Overall, the Gibbs-based sirnulated annealing algorithm 
generates partitions that accurately represent structured 
regions in the image with locally smooth labeling. A 
globally consistent labeling of the partitions was 
sometimes problematic, particularly, for some of the 
smaller regions. However, we initialized the algorithm 
randomly, which is a worst-case scenario. In practice, one 
might consider alternative initialization strategies. Studies 
on additional scenes should be performed. . 
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Multigrid partitioning of hyperspectral imagery using spectral/spatial 

measures of disparity 

Robert S. Rand, Memfeer, and Daniel M, Keenan 

ABSTRACT 

An unsupervised Bayesian approach to partitioning hyperspectral imagery into homogeneous . 
regions is investigated, where spatial consistency is imposed on the spectral content of sites in 
each partition. An energy function is investigated which models disparities in an image that is 
defined with respect to a local neighborhood system using one or certain combinations of the 
spectral angle, Euclidean distance, and/or Kolmogorov-Smimov (classical and mean-adjusted) 
measures. Maximum A Posteriori (MAP) estimates are computed using an algorithm that is 
implemented as a multigrid process to improve global labeling and reduce computational 
intensity. Both constrained and unconstrained multigrid approaches are considered. A locally 
extended neighborhood structure is introduced with the intention of encouraging more accurate 
global labeling. The effort is focused on terrain mapping applications using hyperspectral 
imagery containing narrow bands throughout the 400-2500 nanometer spectral region, and the 
trials of our experiment are conducted on a scene fi-om HYDICE 210-band imagery collected 
over an area that contains a diverse range of terrain features and that is supported with ground 
truth. Quantitative rneasures of local consistency (smoothness) and global labeling, along with 
the class maps, reported in the experiment, denionstrate a clear benefit of this method as 
compared to the commonly used ISODATA clustering algorithm, where the best results are 
achieved with a constrained approach using an energy fiinction consisting of a combined spectral 
angle and Euclidean distance rneasures. 
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Index Terms: Markov Random Field, Gibbs distribution, partitioning, clustering, hyperspectral, 
spectral angle, Kolmogorov Sniirnpv statistic. 

I. Introduction 

The most cornmon unsupervised methods for partitioning multispectral imagery are KMEANS 
and iSODATA [1], [2], where the measure used in determining the cluster membership is the 
Euclidean distance. These clustering algorithms are aiso sometimes used on hyperspectral 
imagery; however, the Euclidean measure, is not particularly responsive to the fine structure, or 
shapes, in high-resolution spectra, and is often overly sensitive to intensity differences. 
Furthermore, these methods do not involve any type of spatial localization or consistency 
operation. A more sophisticated multivariate approach based on the niaximum likelihood 
estimates of first-order and second-order statistics has been developed [3], but relies pii 
dimension reduction methods to operate effectively on the numerous bands inherent with 
hyperspectral data. This method has been tested using matrix factorization for band selection [4] 
but alternative dimension reduction methods include Principal Components Analysis [1] and 
Minimum Noise Fraction [5]. Some related work to [3] investigates a postprocessing method of 
spatial smoothing using the ECHO classifier that is based on a Markov Random Field (MRF) 
model [6]. 

We investigate an unsupervised method of segmenting hyperspectral imagery, where our goal is 
to smoothly label picture elements (pixels) in a scene according to their material type without the 
use of a reference spectral database or the need for supervised training. Our approach will use 
the paradigm of Bayesian image analysis. 
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An extensive naathematical development of the Bayesian framework for Random Processes and 
image analysis, including applications, is presented in [7]. Topics relating to Bayesian image 
analysis are presented in [8], which discusses the Bayesian paradigm, the Gibbs sampler, and 
methods of estimation. Bayesian approaches using Markov Random Field (MRF) models have 
provided algorithms for producing unsupervised partition maps and detecting boundaries on 
imagery containing one band [9]. Gcmm et al [10] use the Gibbs sampler and Metropolis 
methods in their algorithms for the restoration of imagery, partitioning, and boundary detection, 
to compute maximum a posteriori (MAP) estimates. A method to increase the efficiency of 
global optimization within a general MRF framework for image processing problems has been 
developed by Gidas [1 1], the Renormali^ation Group algorithm (RGA), that performs a 
hierarchical, multiscale, coarse-tp-fihe analysis of imagery, for the restoration of degraded 
images. ^ 

In prior work [12], we investigated a spectral mixture process that was conditioned by Gibbs- 
based partitioning. The Gibbsrbased method used an MRF model with a simple energy function 
based on a spectral arigle disparity measure, where MAP estimates were computed using an 
unconstrained multigrid algorithm. The effort showed that significant improvement in spectral 
unmixing results is achieved on a scene containing many fundamental materials, called 
"endmembers," if the scene is partitioned into locally smooth regions, where a smaller set of 
endmembers is attributed to each region. 

In the present investigation, we further develop and test the multigrid Gibbs-based isegmentation 
approach. The effort is focused on terrain mapping applications using hyperspectral imagery 
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containing narrow bands throughout the visible through short-wave infrared spectral region. We 
cpnsidier both constrained and unconstrained multigrid implenientations. We introduce an 
"extended'- neighborhood structure with the intention of encouraging more accurate globial 
labeling without sacrificing local smoothness. We also develop and evaluate a number of 
spectral/ spatial energy functions. The method can be considered an extension, spectrally, of the 
spatial approach used in [9], [10], with additional increased efficiency provided by a multigrid 
implementation that uses a special case of the general course-to-fine grid structure used in the 
RGA approach [11]. The various algorithm implementations are tested using random 
initialization, although, in practice, one can (and probably should) consider some type of 
preprocessing method for initialization; we do so presenting a worst-case initial scenario. We 
use the ISOD ATA algorithni as a control method in our experinient, and also use it as a method 
of initialization in one of the trials of the experiment. The experiment is conducted on a scene 
from HYDICE 210-band imagery collected over an area that contains a diverse range of terrain 
features and that is supported with ground truth. 



II. Mathematical Framework 

A Bayesian model is discussed in this section that provides a framework for the partitioning of 
homogeneous regions in hyperspectral imagery. The approach is model-based. The process 
, called the "partition process," will identify regions that shall be treated as "homogeneous," 

in the sense that each region contains some culturally similar phenomenon (not necessarily a 
single material type). What is actually observed is G, the observed hyperspectral image cube. 
For G there are both spatial and spectral coordinates. 

The approach will be multigrid, and thus we define a collection of lattices that provides a finite 

sequence of coarser and coarser grids: S!p'^ — > x5^^'^ -> ... J^"^ "^"^ \ where 

4"''^{(/CT„ + l,ya„ + l): 0</<^,0$y<^} (1) 

Ng^^ is the number of stages in the sequence, d„ is a parameter called the "grid resolution" 
which determines the spatial sampling of the algorithm and the spatial resolution (coarseness) of 
the labeling process at stage « (/7=.l,...,yV^^) of the multigrid processing, Nj^^^^ is the nuniber 

of rows in an image, and N^^^ is the nuinber of columns, ^p""^ \% the set of sites that will be . 
used in the partition process at a resolution o^. Our grid structure is a special case of the general 
course-to-fine grid structure used in the RGA approach to image processing [1 1], where the 
coarsest level contains all the large-scale features and each finer level contaiiis correspondingly 
smaller-scale features. The term "site" refers to a generic element of a lattice, as specified in (1). 
At ftill-resolution a isite corresponds to a pixel, but at a coarser scale it corresponds to a squa^re 
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block of pixels. The value of o „ is defined to be the number of pixels per side of each block of 
pixels. At the finest resolution =1. We denote the full-resolution image pixel lattice as J), 
i.e., the resolution at which the observations are obtained, where Sj~ Sp^ , 

Let A = {A/ Xi G (400/7/^, 2500 nm\ / = 1, N^^^ denote the full-resolution wavelength lattice, 
where iV^^^ is the nuniber of spectral bands in the hy Hence, for 

example, the fullrresolution coordinates for the observed hyperspectral G are 
Sj % K\s e S^y A e A; therefore, G{s,X) is the observed intensity at site s and wavelength A. . 

The above sets, S^md Sf.''\ have associated graph sfructures which describe, for a given 
element, its neighboring elenients. These specified structures allow us to sequentially formulate 
well-defined probability models for 

Now if 5 is one of S^^"^: S= {s^,. . . ,^^}, where Mis the number of sites on the grid, then an 
associated neighborhood system = {^^, j-g if} is specified (details in Section III A) which is 
a set of lattice graph structures of neighbors, such that s ^ and ssQ^ if and only if r e 
Within this setting, our goal is reduced to finding the value of X'' which maximizes Pr(X'' I G), 
for the observed G. 

is a discrete labeling process which associaites a label with each site, and is treated as an 
MRF. The realizations of vT/ are 4^ e r, where T ={1,2,. . . , A^^a?a} is the set of partition 
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labels. N^^^f^ is the total number of labels assigned to the partitions. The configuration space 
for partitioning is then determined by N^^^^^ and a„ (the label resolution), and denoted by 

The development for JT^ is similar to that of [9] with the important addition of a spectral 
dimension. We treat as a Markov Random Field (MRF) with respect to a graph 
{^S?"?,^^^"^}, where is the entire neighborhood structure for the partition process at 
resolution g^, and is the specific neighborhood at site s for the partition process at resolution 

We exploit the Gibbs equivalence theorem [10] that recognizes the MRF as a random 
process with a Gibbs distribution. The advantage of adopting this approach es that we can now 
compute the maximum a posteriori (MAP) estimate that maximizes the probability Pr(y |^) 
by iteratiyely sampling from the local Gibbs distribution pertaining to each site s e J^'^"^ We 
use 

; pr(< = ^ ix; = ^) = pr(x; =^ i^r e 

. (2) 

' - . ■ ■ ■ .. 

where {X^ = ;r^} means = X and is the 

energy interaction of site s e ^"'^ with the neighborhood (g^. Quite conveniently, the 
computations of are not necessary for purposes of obtaining a MAP estiniate. We denote g to 
be the spatial and spectral realizations of G, which is in contrast to the spatial-only realizations of 
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the simpler G discussed in [10]. Specifically, are the spectral observations at site s, and g are 
the collections of these observations for all the sites. 

The result of this process is a set of partitions (decomposing the set Sf"^): 
5-{s^, qeF, en S'p'^"-}, where E^denptes the collection of all 5 in 5^"^"^ for which ran(iom 
variable Xf takes the value q. For each q, there will be a set, their union forming a 
decomposition of Sp*^'^\ 

III. Description OF AuGORiTHM 

Our proposed algorithm is based on the framework just discussed, We implement a hierarchical 
jEnultigrid approach, where the algorithm typically proceeds in a sequence <7„ = cr^y . .,8,4,2,1. 

We investigate both unconstrained and constrained implementations. 

We shall consider energy fiinctions that respond to spectral disparities in an image using spectral 
angle and/or Euclidean distance, as well as spatial disparities based on the Kolmogprov Smimov 
statistic as in Geman et al [9]. Recall that traditional clustering algorithnis, such as KMEANS 
and ISODATA, only use the Euclidean distance measure to determine cluster membership. 

We noy^ discuss the specifics of the implementation, below, in Sections IIIA-E. We present the 
neighborhood system, define a set of disparity measures and a number of specific energy 
fiinctions of the local Gibbs distribution defined by (2), as well as a inethod for computing the 
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corresponding MAP estimate. We then discuss two niethods for proceeding through the stages 
of the xnultigrid algorithm. Later, in Section IVC, we define perfomiance measures to 
characterize the smoothness of the partition labeling. 

A. Neighborhood System 

The typical MRF neighborhood consists of the closest four or eight neighbors to a site; below, 
we will refer to this as the "compact" neighborhood. However, in order to extend the reach of 
the neighborhood to encourage a more global labeling without a tremendous increase in 
computation, pur algorithm uses an "extended" neighborhood that consists of a predetermined 
near, interrriediate, and far neighborhood structure that still satisfies the (if and only if) 
neighborhood condition, mentioned above. Near, intermediate, and far are relative to a given 
scale (crj. So, the extended neighborhood is a local neighborhood, but "less so" than the 
compact one. 

The structure of the extended neighborhood system fornis a rather simple pattern, The 
neighborhood of a site will consist of near, intermediate, and far neighbors, and they are all at 
specific multiples of Oiir local neighborhood is the eight closest (the surrounding perimeter) 

sites enclosing the site s. The intermediate and far neighbors form a similar pattern as the near 
neighbors but are at a fiirther distance that is a specific multiple of <T^. Figure 1 shows such a 

pattern. The intermediate and far neighbors are used tp aid in obtaining a faster global solution 
and to enable the algorithni to remember label assignments computed by the algorithm at prior 
, stages. Of course, different neighborhood structures will produce, when the conditional 
probability measures (expression (2)) are pieced together, (potentially) different joint probability 
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measures on the overall configuration space. More precisely, the resulting probability measure 
using the "near, intermediate, and far" neighborhood structure is (potentially) a different 
probability measure from that using only near neighbors. Gibbs sampling and/or simulated 
annealing are potentially being applied to two different posterior rneasures. In Section IV, 
experiments are performed, using HYDICE data, comparing the different neighborhoods. The 
"extended" neighborhood appears to produce more accurate global labeling, and does so without 
sacrificing local smoothness. 

B. Disparity Measures 

The energy fiinctions that will be defined momentarily depend on difference measures, called 
"disparity measures j" which indicate the dissimilarity between a site s and a neighbor r. We 
shall investigate disparity measures that measure difference in a number of ways — specifically, 
the spectral angle measure, the Euclidean distance, the classical Kolmogorov Smimov (KS) 
statistic, and a mean-adjusted KS measure. 

For a pair of sites at the site locations r and s, consider a set ^={D^J^\ Dll\ D^^J, 0^^}, where 
the elements of D^ ^ are the spectral angle DjJ], the Euclidean distance D^f,\ the classical KS 
statistic Pll\ and a mean-adjusted KS distance , between the sites s and r. The individual 
distance measures D^'] are defined as follows: 
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Z^^>^sup,^[l//,(^^^^^ (6) 

Equations (5) and (6) need some further discussion. Generally speaking, the KS statistic is used 
to measure the significance of the difference between two univariate empirical distributions. In 
our case, the variable j is the response of a specific band j centered at site r in the 

hyperspectral cube (Note that although we presently just select a single band of imagery, we 
could easily modify the approach such that we use the mean, the maximum, or some other 
consolidated measure of a multiple number of bands), F^grj) deriptes the empirical distribution 
fiinction modeling a window of piixels centered at the site r, and H^(g^ j) denotes an adjusted 

empirical distribution function modeling a window of pixels centered at the site r, where the 
mean of the window of pixels is subtracted from F^(g^j). We denote N^.^ as the number of 

pixels in the spatial window. Therefore, D^^ in (5) is the least upper bound of all pointwise 

differences 1 F^rj)" ^siSsj) Likewise, Pl^^^ in (6) is the least upper bound of all pointwise 

differences \H^{gJ-H^ig^ .)l 
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C. Energy Functions 

For each individual disparity measure in (3) - (6), an energy function is defined as 
Vf/^'U^^^^ i-l4 : (7) 

where the delta function 5(jc^,;c^)=l forx^=jc^ , and (5(jc^,^p==0 for x^^^x^. . 

The ftinction 7)(P^|]) niaps D^'^^ to the set {-1, 1}, where /(D^']) = -l indicates a non-significant 
distance between r and s, and a value /)(D^J^) = 1 indicates a significant distance. If D)^1$:K^jI^^^, 
then f -iPlil) = - 1, and otherwise /.(Dj|,^) = 1 , where the threshold Kjl^,h determines the 
sensitivity of the algorithm to the EJf'l distance measure! The effect of this energy function is to 
accumulate the value of / (Z>^'-) whenever the label of a site s is equal to its neighbor at r. 

For i=]J the value of the threshold Kf^^resh determined through observation; specifically, 
f^rLsh is determined by observing the differences in spectral angle for samples of the sanie 
material verses the differences in spectral angle for samples of different materials, and Kj^r^^h 
determined by observing the differences in Euclidean distance for samples of the same niaterial 
verses the differences in Euclidean distance for samples of different materials. If the imagery is 
calibrated to reflectance, then the value attained should be a fairly universal value, presuming a 
representative set of materials has been analyzed. If the imagery has not been calibrated, then 
this value will depend on the range and scaling of the uncalibrated data. 
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For a probability of false rejection a can conceivably l^e set (say a=.01) and the threshold 
could then be determined by analytical values that characterize the KS distribution (tabulated 
values are available for small windows of A^^^,^ < 45 pixels, and formulas are designated for 
larger windows). Applying the criterion in this fashion corresponds to performing the usual . 
hypothesis testing procedure. As will be discussed further in Section FV, this test of hypothesis 
is too conservative to be useful in practice when applied to the non-adjusted distributions for /-J, 
and experiment observation is again needed to set the t^^^ 

In addition to the energy functions (7) that respond to a single disparity measure, we can 
constmct an energy function that responds to any pair of disparity measures, of any subset 
P^"*" cD^ /of the set of disparity measures J , For example, we could defin^ 
P'f^={Dl}l, D^^^y as the subset that specifies the spectral angle and Euclidean distance between 
sites r and For such subsets, (7) can be generalized as follows: 

: Us{x',g)^'^5(^:\xyrtm (8). 

Appropriate threshold values for calibrated data are given in Section IV, based on observations 
for 1^1,2,3 and based on a statistical test at a=.01 for i^4. 

Conceptually, this energy function is an accumulation of attractive and repulsive interactions 
between sites and neighbors. Whenever a site and one of its neighbors have identical labels they 

are in a common state. In this case, they interact, experiencing an attractive interaction for 

... ^ 

l^r,l- f^Thresh ^^^^ Tepulsivc interactiou for D^,^^>A^^aLa' ^he attractive interaction contributes 



negative energy to (7) or (^), and the repulsive interaction contributes positive energy. 
Whenever a site and its neighbor are not in a comnipn st£^^^^ 

D- CALCUUVTING THE MAP ESTIMATE 

A simulated annealing technique is used to obtain the MAP estimate by the Metropolis algorithm 
as described in [7], [10] and [13]. For a particular site s and iteration we compute the ratio of 
two local Gibbs distributions: , 



r.(>0 = ^*'"''''"^"'^"' (9) ; 



Specifically, Ys(^) the ratio of the Gibbs distribution for a proposed configuration and the 
Gibbs distribution of the current configuration;/', where the in (2) for each distribution 
cancels. The proposed configuration l^is identical to that of , except for a proposal label 
change at site s. The value of Ys(^) "^^^ to determine whether (or not) a transition occurs, in 
which case the configuration of at iteration k+I shall, be set to the proposed one. A transition 
is always made whenever yM)> L If 7/>f) < 1 then the transition is made with a probability of 
y^(^). This latter event occurs whenever yX^) > ^, where ^ is a threshold probability that is 
chosen uniformly in the interval (0,1). Otherwise, no transition occurs (the configuration of at 
iteration will remain unchanged). An annealing schedule is defined through a temperature 
We follow an annealing schedule proposed in [7] and [10], namely, 

/.= ^ (10) 
^ ln(l + >^) K J 



The parameter C corresponds to the niaximum depth of a "well" in the energy fimction (7) or (8), 
and is discussed in [7]. A large value of C starts the computations at a high temperature. A large 
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C is theoretically appealing because a sufficiently high temperature guarantees convergence to 
the niode of tiie Gibbs distribution (a global solution) [10]; however, computationally, a large 
value is undesirable because it results in very slow convergence and long computation times. 
Selecting ai small value of C tends to cool the system too quickly and causes the algorithm to get 
stuck in a local minimum, where convergence to the mode of the distribution is not achieved. 
Because we are using.a multigrid approach, the computing tirne can be reduced considerably by 
using larger values for C at the coarser grids, but smaller values at the finer grids. Practical 
values of C are discussed during the experiments in Section IV. Also as a practical matter, once 
the temperature lowers considerably and a processing stage nears its end, we fast freeze the 
process by dropping the temperature to zero. 

E. Multigrid Approaches 

The multigrid algorithm provides a multi-scale sequence of partitioning with the corresponding 
site lattices (label images) having different resolution cells. The coarse levels contain the 
macroscopic features in a scene, whereas the finest level contains all the microscopic features. 
This is intuitively appealing in the sense that human vision i$ often thought to work this way. 
For example, if a person is looking at a painting from a distance he only sees the larger 
macroscopic feature$. As he proceeds closer to the painting the smaller microscopic details 
became apparent. In order to truly understand the picture, the person may need to view multiple 
scales, and perhaps even step back and forth a few times to "partition" the scene into meaningful 
patterns. 



We investigate two approaches to proceeding through the stages of the multigrid algorithni that 
are consistent with aspects of the Renorrnalization Group Algorithni (RGA) method developed 
by Gidas in [1 1] for the restoration of degraded images. The first approach denoted as 
"unconstrained simulated annealing," uses a grid-coarsening and cascading technique that is 
consistent with one of the techniques used in the RGA niethod. Specifically, the fine-to-coarse 
sequence described in Section II with partition lattices defined by (1) is reversed to a coarse-to- 
fine sequence: Sp"'''^ Sp^ S^pK The results of stage are used to initialize the 
algorithm at the next finer stage However, unlike the RGA method, no constraints are 
applied. 

The second approach denoted as "constrained simulated annealing" uses the same cparse-to-fine 
grid sequence; however, the following constraint is applied: 

. Pr<"'(;c''<'''M.x''<"-')= n W<-. ' . ' 

• ■ ks^! ' ' ■ ■ ^ 

This constraint permanently fixes the labels assigned on completion at stage « as the algorithm 

proceeds to stage n^7. 
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IV Description of Experiment 

A- Data 

A HYDICE scene is used in the experiment that was extracted from data collected over a 
geographic area near Fort Hood, Texas at a 10,0()0 feet altitude, on 16 May 1996. The scene 
consists of 300 samples by 600 lines of pixel vector data with a spatial resolution of 
approximately 1.5 meters. Because of the severe degradations that occur in atmospheric 
windows of the spectrum and sensor degradations occurring in certain bands, a maximum of 1 17 
bands were used to compute the disparities within the energy functions. These 1 1 7 bands of data 
were calibrated to reflectance by the Empirical Line Method based on the known spectra of the 
calibration panels that were located in a region outside this particular scene. 

Ground truth, describing certain ground features in the scene, is based on data collected during a 
site visit to the area during April 2000. The data for this experiment is compiled as a graphic 
oyerlay onto HYDICE Band 49 and is shown in Figure 2, along with associated descriptive 
information listed in Table I. Because of the time difference between scene acquisition and the 
site visit, not all portions of the scene could be documented with certainty. In addition to the 
ground-feature data, spectral data of a diverse range of materials, measured by field 
spectrometers and collected independently of this experiment at other geographic locations, were 
used to determine an appropriate quantization threshold for the spectral angle and Euclidean 
distance disparity metrics, as discussed below. Ground-measured spectral data within the study 
site were not available. 



B, Trial Parameters 

The experiment consists often trials: an ISQDATA run, five unconstrained runs that are 
initialized randonily, one unconstrained run that is initialized by prior ISODATA results, and 
three constrained runs that are initialized randomly. The ISODATA trial uses the iinsupervised 
ISODATA clustering method \yith six classes to be used as a reference against which to test the 
other trials. Table II indicates which of the disparity measures (3) - (6) are used in the energy 
ftmction, where (7) is used as the energy fimction for Trials 1, 3, and 4, and (8) is used as the 
energy function for the remainder of the trials. Trials 2 - 5 are implemented with the parameters 
shown in Table Ilia. These trials proceed in an unconstrained multigrid sequence, with grid 
resolutions a„ = l6, 8,4, 2, 1 . The algorithm is initialized randomly and it then proceeds from a 
grid resolution of = 16 to a, = 1 . As shown in the table, a range C = (1 .0-9.0) was used 
which is below the upperbound discussed in Section HID. The value of C was set to its highest 
value at the initial stage of a run (corresponding to the highest temperature). In addition to the 
practical coniputational advantage, a lower C value allows the algorithm to "remember" the 
results at previous stages. The distances of the intermediate and far neighbors shown in Table 
Ilia, corresponds to fixed multiples of the grid resolution, 4a „ and 8cr^, respectively. 

As indicated by Table II, Trial 6 is initialized by the results df the ISODATA algorithm using the 
same disparity measures and energy function used in Trial 2. 

Trials 7-9 are constrained multigrid runs. These runs are initialized randomly using the same 
disparity measxires and energy function used in Trial 2. As indicated by Table Illb, the 
parameters of Trial 7 are almost identical to tiie parameters in Trial 2. The only differences are 
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that C is higher for some of the intermediate grid sizes along with an appropriately higher 
number of iterations, The constraint should allow us to run the algorithm at a somewhat higher 

temperature (i.e. higher initial C values) without losing the information at coarser grid sizes - we 

■ ■ - ' s ' ■ ■ 

wanted to test this. ' 

Trial 8 is a constrained run that tests the performance of the algorithm using the typical conipact 
neighborhood system, rather than the extended neighborhood system we propose aiid have tested 
in the prior runs. As shown in TableTIIc, Trial 8 is identical to Trial 7 except for the compact 
neighborhood system. Again, by compact, here we mean the eight nearest neighbors. 

Trial 9 is a constrained run that also tests the performance of the algorithm using the typical 
compact neighborhood system. Because of the anticipated loss in global information using this 
neighborhood system, we attempt to compensate for the compact neighborhood system by 
starting the algorithm at a coarser grid size than that used in Trial 8 (or any of the prior trials). 
Instead of the five stages used in the other trials, this trial is run with six stages using the 
parameters shown in Table Illd. 

C. Measures OF Performance 

Quantitative results are reported for the trials, making specific use of the ground truth (GT) 
infonnation to generate confusion matrices and to compute a number of quantitative measures of 
performance. The quantitative measures include global indicators of labeling accuracy such as 
commission errors (percentage of class labels assigned to a partition that correspond to 
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categories other than the one attributed to that partition) and omission errors (percentage of GT 
labels that are assigned to partitions not associated with a specific GT category), and indicators 
of labeling smoothness ^uch as clutter and speckle to indicate labeling smoothness. 

CI Global Measures 

Whereas, errors of omission and conimission are widely used to report the performance of 
supervised classification algorithms and are rather straightforward to compute for these 
algorithins, the computations of such errors for unsupervised classification are not quite as 
straightforward. In supervised classification, the identity of the classification labels with respect 
to ground truth is typically determined in advance and so there is an exact correspondence 
between the identity of the labels and the ground truth that is built into the confusion matrices. 
However, in unsupervised classification, the clusters develop independently of the ground truth 
and the correspondence between the identity of the classification labels and GT labels must be 
made afterward, and is often not quite so straightforward. lii this experiment, a correspondence 
of the partitioning labels with three GT categories (Trees, Grass, and Asphalt) was attempted. A 
successfiil correspondence was usually (but not always) made for assigning three out of the six 
partitioning labels to these three categorieis. - 

C2. Local Measures To Characterize Classification Smoothness 

Smoothness in classification labeling can be measured from the viewpoint of clutter and speckle. 
Suppose we have a region in a scene that is known to be homogenous. Therefore, all the pixels 
in this region should be assigned to one label. Clutter can be defined as the number of pixels 



assigned to labels that are different from the predominant label in the region. In addition, if there 
are a significant number of pixels with differing labels in a region, these differing labels could be 
scattered in a speckle-type pattern with a rough texture or they could be clumped together 
.exhibiting smooth texture. In order to characterize the roughness of the labeling texture in a 
region we will define measures of what we call *'specW^^^ 

One clutter-type measure to characterize smoothness (or homogeneity) of classification is easily 
computed from the typical confusipn matrices generated to measure classification performance. 
The clutter measure, Y^autter d^fii^^d as 

J^^'^^= At-I'^^^^/^r} ' /VI.A'^ (12), 

where Y^^ is the number of sites assigned to the label in the f ground truth (GT) region, 

■> 

^GTreg nuiuber of GT regions being evaluated and N^j.^.^^ = ^ ^cTsUes - The quantity 
yV^^^^ is the number of sites contained in the Z'* GT region and /= 1, A^cn^- 

As an indication of the roughness of the labeling in a homogenous region, we consider a measure 
of speckle that counts the number of neighborhood labels that are different from each site in a 
region. The measure K^^^]/^ is defined with respect to a* very simple neighborhood 
consisting of two neighbors - the neighbor to the left and the neighbor above - of a site: 
K^.= I Z (13) 

where *y^^ is the set of ground truth sites. 
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A second speckle measure ^sp^^k ^^^^ ^® defined with respect to a more symmetric 
neighborhood ^^"""^ • consisting of four neighbors - the neighbors to the left and right, as well as 
the neighbors above and below ^ of a site: 

K<^-I L (l-^J (14) 

V. Results 

Figure 3 and Table IV show the results of the ISOD ATA clustering method applied to the image 
data for six classes. The comrnission errors for Trees, Grass, arid Asphalt are 22.48%, L32%, 
and 68.66%, respectively. Although ISODATA provides a fairly good partitioning of the scene, 
it is overly sensitive to minor vegetative differences in the forested areas, denoted by "D1-D6'' in 
the ground-truth overlay (Figure 2) as well as by shade. This is apparent with the high 
commission errors for Trees and Asphalt. Omission eirors also indicate the same problem. In 
particular, omission errors show tha^t another label associated with Al (Asphalt) is prevalent 
throughout subregions of DI -D6 that can be attributed to shade. Specifically > 27% (1015 out of 
3764) of the P1-D6 pixels are labeled with the sanie cluster associated with AsphaU. Also 
15.7% (559 out of 3559) of the GNG4 pixels are labeled with the same cluster as associated with 
Trees. Furthermore, the clustering does not provide a smooth and spatially consistent mapping 
of regions, and exhibits ^ significant amount of speckling. Measures of smoothness exhibited the 
values K^/^,,^^ = 0.22, K^^^^]/^ = 1641, and K^^^^j^^ = 3359. These measures will be compared to 
those of the Gibbs runs at the end of this section (Tables XIII and XIV). 
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Using observations of spectxal angle distances between the vatious .pairs of materials we set 
^Thresh = 1 1.0, which remained constant throughout those partition trials involving the angle 
measure. ^Using observations of Euclidean distances between the various pairs of materials, we 
set K^^^^^^f, = 100.0, which also remained constant throughout those partition trials involving that 
measure. For both measures and most of the pair-wise signature combinations, there is an 
obvious gap between similar materials and different materials, except for combinations involving 
a vegetation signature. Unfortunately, vegetation has a fairly large variance, which makes the 
threshold boundary a bit fuzzy. During these trials, we observed that the exact value of these 
thresholds did not seem to be very critical. 

Our initial thought was that the values of K^T^^^^h ^z^LA^^^ld be selected analytically based 
on the theoretical properties of the KS statistic. However setting the value of K^rhresk tP any of the 
commonly used a values (probabilities of false rejection) resulted in the algorithm finding a 

significant difference between alniost everything. Consequently, the algorithm using the 
classical KS statistic based on standard a values could not generate a partition map containing 

any of the structure in the scene. Our inteipretation is that the KS statistic is very conservative 
about accepting a hypothesis that two distributions are similar, and that most of the observed 
intiensity differences were statistically significant. After a number empirical runs, we observed 
that meaningful structures only occurred at KS distances well beyond the usual theoretical 
Values, and we determined empirically that a value of AT^^^^ = 0.5 was suitable for a spatial 
window size of 5x5. 
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Figure 4 and Table V shows the results of Trial 1 testing the spectral angle measure in the energy 
function! This trial was successful in characterizing the structure of many of the natural and 
cultural features in the scene without much of the undesirable clutter prevalent in the ISODATA 
scene shown in Figure 3. For example, the algorithm was not overly sensitive to minor shade 
differences such as in the forested areas, denoted by "D1-D6" in the ground-truth map. 
However, the algorithm was also insensitive to other differences that we would like to 
discriminate. For example, there is missing structure within the road mtersection, denoted by 
"Vl," and the rooftops of the apartment buildings (asphalt shingles) denoted by "J2" are 
confused with the surrounding asphalt parking lot "A3." The, commission errors for Trees, 
Grass, and Asphalt are 4.9%, 2. 1 6%, and 33,07%, respectively. Measures of smoothness 
improved over the ISODATA trial, exhibiting the values K^,,,,,, = 0.017 , K^^^^4 = ^22, and 

^Speckle . 

Figure 5 and Table VI shows the results of Trial 2 using both the spectral angle and Euclidean 
distance measures in the energy function. A coniparison of Figures 4 and Figure 5 shows a clear 
advantage to using the two measures in the energy function. Many of the cultural structures 
missed in Trial 1 are detected in this trial without introducing undesirable vegetation clutter. For 
example, the structure in "Vl" is now evident and so are the rooftops in "12." The commission 
errors for Trees, Grass, and Asphalt are 24.08%, 3.88%, and 20.40%, respectively. There has 
been a significant drop in the commission error of the partition associated with Asphalt 
compared to the ISODATA trial. Regarding omission error. Table VI shows that the regions Dl- 

D6 are labeled quite homogeneously and are not adversely affected by shade, as was the case 

^ ■ ■ ' . ■ ' 
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with ISODATA. Specifically, 0.9% (34 out of 3764) of the D1-D6 pixels are labeled with the 
same cluster associated with Asphalt. Unfortunately, there is a significant problem with the 
grass region G4, where 8 1 8 out of the 832 GT sites were assigned the label associated with 
Trees. Measures of smoothness have improved over the ISODATA trial, with values 
. K^^^^ 0.03, K^^.=. 395, and K<^^ 795. 

Figures 6 and 7 show the results of the single-measure KS trials using a 5x5 spatial window. 
Figure 6 shows the results of Trial 3, where the classical KS statistic is tested using an 
empirically-derived value K^fl^^^^ = 0.5, as mentioned above. Trial 3 was the best performing 
trial over a number of other similar experimental trials at 5x5 windows (not shown) using 
different values of K^^lg^h- Nevertheless, the results poorly represent the structure in the scene. 
Figure 7 shows the results of Trial 4, where the mean-adjusted KS statistic is tested using an 
analytically-derived value of K^H^^f^, as mentioned above. Trial 4 provides a very reasonable 
partitionitig of the scene, particularly when considering that only one band of imagery is used in 
the analysis. This measure detects some structure in the building denoted by feature "B l" in the 
ground-truth map, not found by the other measures, as well as the cars in the parking lot of 
feature "Al." However, the mean-adjusted KS measure missed a number of cultural features and 
was unable to distinguish between Grass and Trees. Edge effects are also prevalent throughout 
the partition map. Not shown are the results of similar trials conducted using a 7x7 spatial 
window that generally produced inferior results compared to the 5x5 window. 

Figure 8 and Table VII shows the results of Trial 5, which combines the spectral angle, 
Euclidean distance, and mean-adjusted KS statistic measures in the energy function. The 
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partitioning exhibits more detail than any of the other measures; however, edge effects are 
prevalent and give the impression of too much clutter, particularly in the cultural regions. 
Measures of smoothness have degraded compared to the previous Gibbs trials (Trial 1 and Trial 
2), exhibiting values of Ka,„er = 0-07 , K^^;^^,, = 665 , and K^^l^ = 1377. This measure might 
prove to be usefiil; in the future, if the energy fimction could be modified to discourage the 
development of edges. The use of penalty tenris in the energy function, perhaps of the form 
suggested by Geman et, al. [9] might be beneficial in removing such artifacts. On the other 
hand, the tendency for this measure to produce edge artifacts might be iii fact be suggesting an 
alternative use for this energy function as a measure in an edge inducing algorithm — perhaps a 
Gibbs-based boundary detection algorithm. 

All the above trials for Figures 4-8 were performed using identical control parameters, as listed 
in Table Ilia, in order to provide some consistency for comparing the different methods. The 
results indicate that an energy function using the combined spectral angle arid Euclidean distance 
measures is the preferred method, unless an energy function involving the KS measure could 
somehow be modified to discourage the development of unwanted edges. Based on these results, 
the subsequent trials (Trials 6-9) were focused on a specific energy function - the one that 
combines the spectral angle and Euclidean distance measures. 

Figure 9 shows the partition map resulting from Trial 6, where the Gibbs algorithm is initialized 
by the results from the ISODATA algorithm, and Table VIII lists quantitative resuhs. The 
commission errors for Trees, Grass, and Asphalt are 6.09%, 3.43%, and 20.19%, respectively. 
So there has been a significant drop in the commission error of the partition associated with 
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Asphalt, as well as for Trees, compared to the ISODATA trial. Measures of smoothness have 
improved compared to the ISODATA trial, with values of Kc,„„,, - 0.0334 , K^p^^^^ = 372, and 
^^speJt/e ^ 794. These results show that Trial 6 provides substantial improvement over the results 
of the standard ISODATA trial (Figure 3 and Table IV). 

We now shift attention to studying some additional properties of the algorithm. Given that the 
combined spectral angle and Euclidean distance measure has shown the most promise, this 
aspect is fixed for the remaining trials, and we now look at constraining the multigrid process 
and also the examine the effect of using the naore typical compact neighborhood structure in the 
energy function. 

iFigure 10 and Table IX show the results of Trial 7. As mentioned in Section I VB, this trial is 
similar to Trial 2 except that it that it is constrained and it is run at higher initial temperatures 
during the intermediate stages. A comparison of the results with Figure 5 and Table VI show a 
better global labeling using the constrained algorithm. The commission errors for Trees, Grass, 
and Asphalt are .94%, 3.03%, and 21 .1%, respectively. The omission errors also improved 
compared to Trial 2. Specifically, 1 .5% (57 out of 3764) of the P1-D6 pixels, are labeled with 
the same cluster associated with Asphalt. Also 0.1% (4 out of 3559) of the G1-G4 pixels are 
labeled with the some other cluster and the large omission error for G4 in Trial 2 has been 
reduced to 1 .7%. Measures of snioothness for Trial 7 were slightly degraded compared to Trial 
2, exhibiting values of Kc^„„ = 0.0424, K<,^,^, ^ 
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Figure 1 1 and Table X show the results of Trial 8. This tri^ 

parameters as Trial 7, except that the typical cornpact neighborhood is used in the energy 
function rather than our proposed one. A comparison of these results with Figure 10 and Table 
IX shows a significant degradation in global labeling. For the trees, in Trial 8, the GT regions 
Dl ; D3, and P4 are assigned different labels than D5-D6. For the grasses, in Trial 8, the GT 
regions for G1/G3, G2, G4, G5, and G6 are assigned different class labels. The commission 
errors for Trees, Grass, and Asphalt are 21.45%, 34,08%, and 81.32%, respectively. Measures of 
smoothness for Trial 8 were slightly improved oyer Trial 7, exhibiting values of l^ciuner = 0.0401, 
ICS-467,andK^= 930. ^ ■ 

Figure 12 and Table Xl show the results for Trial 9. This trial was made using the same 
parameters as Trial 8, except for the addition of one coarser stage in tiie multigrid process. The 
addition of this stage was made in an attempt to extend the reach of the compact 5x5 
neighborhopd by beginning the hierarchical processing at a coarser scale. Visually, the class 
map seems slightly better than the one for Trial K. However, the global labeling is still 
significantly degraded from the results achieved in Trial 7 or Trial 2. In fact; the task of making 
a correspondence between a GT label and a partition label for purposes of establishing 
commission and omission errors becomes a problem for Grass and Asphalt. Indeed, for purposes 
of commission error we need to treat LBL03 as both Grass and Asphalt, which is troublesome 
and makes the usefulness of this measure questionable for this trial. The cornmission errors for 
Trees, Grass, and Asphalt are 39.03%, 27.82%, and 80.26%, respectively. Measures of 
smoothness for Trial 9 were almost identical to Trial 8, exhibiting values of Kauner - 0.0483, 
4^i= 477,andK<-l= 924. 
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Table XII provides a summary of commission errors for three of the predominant GT categories 
and for eight of the trials. This table shows that the top two trials in terms of global labeling 
performance are Trials 6 and 7, these trials both used the combined spectral angle and 
Eyclidean distance measures in the energy function. The global labeling of Trial 6 was helped 
by a non-random initialization. Trial 7 was Helped by the use of the constrained multigrid 
implementation. 

Tables Xin and XIV summarize the homogeneity measures. There is clearly an improvement in 
homogeneity over the ISODATA algorithm for the Gibbs-based trials using the proposed 
neighborhood structure shown iri Figure 1 (Trials U 2^ 5, 6, and 7). The trial using a compact 
5x5 neighborhood also shows notable improvement in homogeneous labeling (Trial 8). 

The far right column of Table XIV lists the ratio of K^^^^J^}^ and ^^specL a means to examine 
any correspondence between the two measures. This column show's K^^^^^)^ and K^^^jti ^F^ 
different by a scale factor of approximately two, which is not surprising. 

Finally, the results of Trial 7 (the best performing trial) are shown in Figure 14 as a color class 
map, with Figure 13 showing a red, green, blue (RGB) true-color composite of the scene. 
Visually, the color map confirms the results discussed above, showing excellent globally-smooth 
classification labeling for the majority of the regions. The difficulty that the method sometimes 
has with narrow and small features is also apparent. 
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V Conclusions 

A number of energy functions were investigated involving the use of the spectral angle measure, 
the Euclidean distance, the classical Kolmogorov Smimov (KS) statistic, the mean-adjusted KS 
measure, and certain combinations of these measiires . The results of this experiment indicate 
that the energy function using a combination of the spectral angle and EucUdean distance 
nieasures will produce excellent results. This energy function is the preferred one in the present 
context. ' 

The tendency for the KS measure (as a single measure, or in combination with other measures) 
to produce edge artifacts may actually indicate an alternative use for this energy function as a 
rneasure suitable for an edge-inducing algorithm — perhaps a Gibbs^based boundary detection 
algorithm. It should be reiterated that the KS measure was calculated at a single wavelength. A 
current avenue being explored is where the KS statistic is calculated as a voxel neighborhood 
(two spatial and one spectral). 

Overall, the Gibbs-based algorithm generates partitions that accurately represent structured 
regions in the image with locally smooth labeling. For most of the bigger regions, pixels with 
the same label represent the same type of phenomenon on the ground, globally across the image, 
and regions with different labels correspond to different phenomenon. A globally consistent 
labeling of the partitions was sometimes problematic, particularly for some of the sm^U as well 
as narrow regions. This problem can be reduced significantly by using (1) non-random 
initialization and/or (2) the constrained muhigrid approach. This is explicitly demonstrated in 
Trials 6 and 7, respectively. The quantitative measures of local consistency (smoothness) and 
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global labeling, along with the class maps demonstrate a clear benefit of the Gibbs-based method 
as compared to the traditional ISODATA clustering algorithm. 

The proposed neighborhood structure discussed in IIIA was shown to provide significantly better 
global labeling than the typical MRF compact neighborhood without any significant degradation 
in labeling smoothness. A comparison of the results of Trials 8 and 9 with the corresponding 
trials using our proposed non-compact neighborhood structure explicitly niake this point. 

. . ^ . ^ .......... ^ . . 

The Gibbs-based algorithm can be used as a stand-alone method; however, it is also envisioned 

to work within an overall framework involving other processes. In particular, the partitioning of 

hyperspectral imagery is viewed as a mechanism to condition a spectral mixing process as 

proposed in [12]. Under this framework, the ultimate identification of materials (and the 

corresponding abundances) in a sceiie is a consequence of the spectral mixing process, so that a 

precise global labeling is not required at the partitioning stage. 
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Table 1. Description ofthe ground truth (GT) information, as displayed in Figure 2. 



GT Label 


Descriotion of the Area 


Al 


Asphalt parking lot, newly paved.. 


A2-A3 


Asphalt parking lot. A3a and A3b opposite sides of apartment complex, different ages. • 


A4-A5 


Asphalt road 


Bl 


Rooftop of a department store, light-gray asphalt shingles, skylights, yellow gas pipes. . . 


B2a,b 


School, rooftop mostly coated bubbly light- gray material (B2b), top edge is gravel (B2a). 


CI 


Concrete pavement to loading dock, some tire marks. 


C2a,b 


Concrete parking lot surrounding small bank building. C2b is brighter than C2a. 


bl 


Trees: Black Willow, Texas Sugar Berry, Dogwood, Texas oak, Elbpn bush, Green Bryer. 


D3 


Trees 




Trees: Texas Oak, Burmelia, and Texas Hackberry. 


" tSc 
UD 


irees. Kea rsuQ, lexas uaK, ana Kea puo. 


D6a,b,c 


trees - Deciduous holly, Bois D'Arc, Osage Orange, Texas Oak, and Bow Wood 


"D7 


Scattered junipers throughout the grassy areas. 


El 


Church and school with metal rooftop. 


E2-E6 


Corrugated steel buildings. 


G1-G6 


Tall wild grass. G5 and G6 near drainage and healthier. 


HI 


Healthy maintained grass near bank building. 


Jl 


Private residence, stand-alone! 


J2 


Apartnient buildings Witli asphalt-shingled rooftops. 


J3 


Church with dark brown asphalt-shingled rooftop. 


LI 


Shade of Apartment building. 


Ml 


Tennis court, concrete with blue/gray coating, and lots of dirt. 


". Ni 


Playground area with exposed soil. 


01 


Encircled area is residential with moderately-priced homes and concrete driveways. 


Rl 


Light linear pattern is exposed bedrock. 


R2-R3 


Exposed soil. 


VI 


Asphalt road intersection with exposed soil along the shoulders. 


Wl 


Playground - grass with exposed soil, gravel, and some asphaU. 
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Table II. Descriptive summjEiry of the trials. 



Trial 


Disparity nieasure(s) 


Description 


1 


Spectral angle ■ 


Unconstrained, 5 stages , 


/2 ■ 


Spectral angle, Euclidean distance 


Unconstrained, 5 stages 


3 


Classical KS statistic (5x5 window}^ 


Unconstrained, 5 stages 


4 


Mean-adjusted KS statistic (5x5) 


Unconstrained, 5 stages 


"J ^- ■ 


Spectral angle, Euclidean distance, mean-adjusted KS statistic 


Unconstrained, 5 stages 


6 


Spectral angle, Euclidean distance 


Initialized by ISODATA 


7 


Spectral angle, Euclidean distance 


Constrained, 5 stages 


8 


Spectral angle, Euclidean distance (5X5 compact 
neighborhood) 


Constrained, 6 stages 


9 


Spectral angle, Euclidean distanceX5X5 CO 
neighborhood) 


Constrained, 6 stages 
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Table III. Control parameters for the trials. The sequence of numbers in the "Dist of neighbors" column corresponds to 
the distance from the center pixel (the "site") to the perimeter, for the surrounding liear neighbors, the intermediate neighbors, 
and the far neighbors. The sequence of numbers in the "# of neighbors" column corresponds to the number of near neighbors,' 
intermediate neighbors, and far neighbors in the neighborhood system. 



Table Ilia. Control parameters for Trials 1, % 3, 4, 5 



Q 


C 


# iterations 


Dist. of nei&hbors 


# of neifihbors 


16 


9.0 


50,000 ^ . ' 


1 6 64 1 28 




8 


3,0 


2,000 


8,32,64 


8, 8, 8 


4 


3.0 


. 2,000. 




8, 8, 8 


2 


2,0 


700 


2,8,16 


8, 8, 8 - ■■ 


1 


1,0 


350 


1,4,8 


8,8,8 


Table Illb. Control parameters for Trial 7. 






G 


c 


# iterations 


Dist. of neiehbors 


# of neighbors 


16 


9.0 


50,000 


1 6 64 1 28 


8, 8, 8 


8 


6.0 


10,000 


8,32,64 


8, 8, 8 


4 • 


3.0 ■ 


2,000 


4 J 6,32 


8, 8, 8 


2 


3.0 


2,000 


2,8,16 


8, 8, 8 


1 


1.0 


350 


1 ft) o 


8 8 8 


Table IIIc. Control parameters for Trial 8. 








c 


# iterations 


Dist. of neighbors 


# of neighbors 


16 


9,0 


50,000 


compact 


24 - . . 


8 


6.0 


10,000 


compact 


24 . 


4 


3.0 


2,000 


compact 


24 




3.0 


2,000 


compact 


24 


1 


LO^ 


350 


compact 


24 . 


Table Illd. Control parameters for Trial 9. 






0 


c 


# iterations 


Dist. of neighbors 


# of neighbors 


32 


9.0 


50,000 


compact 


24 


16 


9.0 


50,000 


compact 


24 ■ ; :; 


8 


6.0 


10,000 


compact 


24 


4 


3.0 


2,000 


conipact 


24 


2 


3.0 


2,000 


compact 


24 


1 


1.0 


350 


compact 


.24 
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Table IV. Confusion matrix off results for the ISODATA algorithm. Labels in the rows are the GT labels, and 
labels in the columns are the partitioned labels. As an example of how to read the table, the region of ground truth 
extracted for CI consisted of 120 pixels. Of the 120 pixels, 13 pixels were assigned to cluster label ISOr04, and 107 pixels were 
assigned to cluster label ISO-05. 



Labels 


ISO-1 


lSO-02 


ISO-03 


ISO-04 


iSO-05 


ISO.06 


NUL^. 


Al . 


0 


0 


, 518 


0 


0 


0 




Bl 


0 


0 


0 


278 


0 


0 


0 


CI 


0 


0 


0 


13 


107 


Q 


0 


C2a 


0 


0 


0 


'' 0 


13 


0 


0 


C2b 


0 


0 


0 


Q 


11 


0 


, 0 


B2a 


0 


0 


0 


• 2 


15 


0 


0 


B2b 


0 


0 


0 


12 


110 


. 0 


0 


Ml 


0 


0 


0 


43 


0 


0 


0 


J3 


0 


0 


33 


44 


0 


0 


0 


A2 


0 


.4 


0 


51 


0 




0 


A3a 


0 


0 


0 


16 


0 . 


0 


0 


A3b 


0 


0 


0 


12 


0 


0 


0 


El 


0 


0 


0 


0 


. 29 


0 


0 


E3 ■ 


0 


0 


0 


11 


13 


0 


0 


E4 


0 


0 


0 


0 


15 


0 


0 


E5 


0 


0 


0 


0 


71 


0 


0 


E6 


0 


0 


0. 


22 


8 


0 


0 


J2 


0 


0 


84 


0 


0 


0 


0 


Nl 


0 


0 


0 


5 


88 


0 


0 


R3 


0 


0 


0 


0 . 


56 


0 


0 


Dl 


203 


0 


68 


0 


0 


1 


0 


03 


231 


8 


138 


0 


0 


9 ^ 


0 


D4 


859 


11 


377 


0 


0 


98 ^ 


0 


D5 


262 


6 


108 


0 


0 


54 


0 


D6a 


279 


4 


43 


0 


0 


233 




D6b 


346 


5 


261 


0 


Q 


41 


0 


D6g 


96 


2 


20 


0 


0 


1 


0 


Gl 


Q 


504 


0 


0 


Q 


8 


0 


G2 


87 


969 


0 


0 


0 


0 


0 


G3 


13 


1143 


3 


0 


0 


0 


Q 


HI 


0 


0 


0 


0 


0 


33 


0 


Wl 


0 


1 


0 


65 


0 


0 


0 


G4 


459 


373 


0 


0 


0 


0 


0 


G5 


101 


0 


0 


0 


0 


71 


0 


G6 


0 


10 


0 


0 


0 


176 


0 


Total GT 


2936 


' 3040 


1653 


574 


536 


725 


0 


Total non-GT 


27775 


55909 


29245 


26055 


10402 


21150 


0 
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Table V. Confusion matrix of results for Trial 1 - the spectral angle disparity measure. Labels in the rows are 
the GT lab.els, and labels iii the columns are the partitioned labels. 



Labels 


LBL-01 


LBL-02 


LBL-03 


LBL-04 


LBL-05 


LBL-06 


NULL 


Al 


0 


0 


518 


0 


0 


0 


0 


Bl 


0 


0 


0 


277 


1 


0 


0 


CI 


0 


0 


120 


0 


0 


0 


0 


C2a 


0 


0 


7 


P 


0 


0 


0 


C2b 


0 


0 


11 


0 


0 


0 


0 


B2a 


0 


0 


0 


17 


0 


= 0 


0 


B2b 


0 


0 


0 


122 


0 


0 


0 


Ml 


0 


0 


0 


2 


41 


0 


0 


j3 


0 


0 


0 


0 


0 


77 


0 


A2 


0 


0 


55 


0 


0 


0 


Q 


A3a 


0 


0 


16 


0 


0 


0 


0 


A3b 


0 


0 


12 


0 


0 


0 


0 


El 


0 


0 


0 


29 


0 


0 


0 


E3 


0. 


0 


0 


0 


0 


24 , 


0 


E4 


0 


0 


. 0 


15 


0 


0 


0 


E5 


0 


0 


0 


0 


71 


0 


0 


£6 


.0 


0 


0 


0 


2 


28 


0 


J2 


0 


0 


84 


0 


0 


0 


0 


Nl 


0 


1 


0 


92 


0 


0 


0 


R3 


0 


0 


56 
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0 


0 


0 


Dl 
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0 


0 


0 


0 


0 


0 


D3 
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27 
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6 


3 


0 


U4 
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2o 


/ 


0 
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Z 


0 


D5 


412 


6 


1 


6 


0 


5 


0 


D6a 


551 


3 


2 


2 


0 


1 


0 


D6b 


620 


10 


7 


10 


1 


5 


0 


P6c 


117 


0 


2 


0 


0 


0 


0 


Gl 


0 


0 


0 


512 


0 


0 


0 


G2 


0 
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0 


1 


0 


1 


0 


G3 


0 
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0 


1 


0 


0 


0 


HI 


0 


33 


0 


0 


0 


0 


0 


Wl 


0 


66 


0 


0 


0 


0 


0 


G4 


15 


817 


0 


0 


0 


0 


0 


G5 


172 


0 


0 


0 


0 


0 


0 


G6 


0 


185 


. 0 


0 


I 


0 


0 


Total GT 


3815 


3386 


898 


1088 


125 


146 


0 


Total non-GT 


41033 


65311 


30630 


13467 


9636 


.5009 


5456 
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Table VI. Confusion matrix of results for Trial 2 - the combined spectral angle/Euclidean distance disparity 

■ \ .• ' - . . 

.measure. Labels in the rows are the GT labels, and labels in the columns are the partitioned labels. 











I Rf 0/1 


r nr a< 


1 R1 HA 




A 1 
Al 


0 


0 


518 


0 


0 


0 


0 


ni 

Dl 


0 


0 


0 


278 


0 


0 


0 




0 


0 




0 


120 


.0 


0 


C^a 


6 


0 


0 


0 


0 


1 


6 




11 


0 


0 


0 


0 


^ 0 


0 


pza 


.0 


0 


0 


0 


17 


0 


0 




0 


0 




0 


122 


0 


0 


iVll \ 


0 


0 


41 


0 


2 


0 


• 0 


J3 


0 


0 


77 


0 


0 


0 


0 


AZ 


0 


0 


55 ' 


0 


' P 


0 


0 


A3a 


0 


0 


16 


0 


0 


0 




AJD 


^ 


p 




0 


0 


0 


p 




9 




0 




19 


0 


0 


VI. 


0 


Q 


0 


7 


0 


17 


0 


11.4 


0 


15 


0 


0 


0 


p 


0 


E5 


29 


0 


0 


0 


0 


42 


0 




2 


0 


0 


24 


0 


4 


0 






p 


0 


0 


84 


0 


0 


Nl 


0 


I 


0 


0 


92 


0 


0 


R3 


56 


0 


0 


0 


0 


0 


0 


Lli 


272 


0 


6 


0 


0 


6 


0 




345 


25 


0 


4 


12 


0 


0 


D4 


1278 


32 


8 


14 


6 


7 


0 


D5 


394 


15 


13 


4 


2 


2 


0 


D6a 


527 


14 


5 


4 


4 


5 


0 


D6b 


601 


10 


8 


12 


11 


11 


0 


D6c 


117 


0 


0 


0 


0 


2 


0 


Gl 


0 


512 


0 


• Q 


0 


0 


0 


G2 


0 


1054 


0 


1 


0 


1 


0 


G3 


0 


1159 


0 


0 


0 


0 


0 


HI 


0 


33 


0 


0 


0 


- 0 


0 


Wl 


0 


66 


0 


0 


0 


0 


0 


G4 


' 818 


0. 


,0 


0 


14 


0 


0 


G5 


0 


0 


0 


0 


172 


0 


0 


G6 


185 


1 


0 


0 


0 


0 


0 


Total GT 


4655 


2939 


755 


352 


680 


98 


6 


Total non-GT 


41394 


64214 


20247 


15298 


14026 


9904 


5453 
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Table VII. Confusion matrix of results for Trial 5 t the combined spectral angle, Euclidean distance, and 
mean-adjusted KS disparity measures. Labels in the rows are the GT labels, and labels in the columns are the 
partitioned labels. 
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V 


n 
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n 
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z 


A 


1 
I 


ft 


Dl 


■ 0 


0 


0 


0 


0 


272 


0 


D3 


344 


25 


0 


4 


5 


8 


0 


D4 


' 20 


9 


1275 


10 


10 


21 


0 


P5 


9 


11 


4 


2 


8 


396 


0 


D6a 


521 


4 


12 


1 


5 


16 


0 


D6b 


597 


5 


17 


10 


.10 


14 


0 


D6c 


115 


0 


0 


2 


1 


1 


0 


Gl 


1 


12 


494 


0 


4 


1 


. 0 


G2 




1046 


1 


6 


0 


2 


0 


G3 


61 


0 


1097 


0 


0 


1 


0 


HI 


1 


0 


32 


0 


0 


0 


0 


Wl 


0 


0 


66 


0 


0 


0 


0 


G4 


0 


824 


8 


0 


0 


0 


0 


G5 


0 


1 


166 


0 


p 


5 


0 


G6 


4 


83 


0 


1 


0 


98* 


0 


Total GT 


1821 


2237 


3737 


350 


251 


1062 


0 


Total non-GT 


26702 


37175 


38199 


19193 


22217 


.21544 


5512 



43 



Table VIII. Confusion matrix of results for Trial 6 - the combined spectral angle and Euclidean distance 
disparity measures, when initialized by the ISODATA algorithm. Labels in the rows are the GT labels, and 
labels in the columns are the partitioned labels. 



Labels 


T Df fkt 

LBL-lll 


L.I5L.-U2 




¥ Dt A>4 
LfDL#-U4 


¥ DT AC 


¥ D¥ A/: 


TVT¥T¥ T 


Al 


0 


0 
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0 


0 


0 


A 

0 


Bl 


0 


0 


p 
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0 


0 


0 


CI 


0 


p 


0 


0 
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0 


p 
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5 


2 


0 


P 
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0 


cib 


1 1 


p 


0 


0 


P 


0 


0 
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0 


0 


P 


0 


17 
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0 


A 

0 


B2b 


0 


0 


0 


1 
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A 


0 


Ml 


0 


0 


P 
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0 
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J3 
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0 
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ri 
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A2 
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0 
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A 

0 
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■A 

I) 
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29 
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0 


P 
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25 


0 


J2 


0 


0 


84 


0 


0 


0 


0 


Ml 


u 


i 


U 


A 


yz 


A 
U 


u 


R3 


0 


0 


A 

0 


r\ 
U 


56 


A 


A 

u 


U I 




0 


Q 


0 


0 


0 


0 


D3 


, 343 


. 27 


I 


0 


I 


14 


0 


P4 


1284 


22 


15 


7 


4 


13 


0 


D5 


392 


14 


3 


.7. 


9 


5 


0 


D6a 


524 


6 


11 


6 


5 


7 


0 


D6b 


598 


9 


7 


12 


8 


19 


0 


D6c 


117 


0 


0 


0 


2 


0 


0 


Gi 


0 


512 


0 


0 


0 


0 


0 


G2 


0 


1054 


0 


0 


0 


2 


0 


G3 


0 


1159 


0 


0 


0 


0 


0 


HI 


p 


0 


0 


0 


33 


0 


0 


Wl 


0 


66 


0 


0 


0 


P 


0 


G4 


22 
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0 


0 


0 


0 


0 


G5 


172 


0 


0 


0 




0 


0 


G6 


0 


3 


0 


0 


0 


183 


0 


total GT 


3759 


3732 " 


718 


466 


477 


312 


0 


Total non-GT 


36442 


71392 


18126 


15793 


13671 


15112 


0 
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Table JX. Confusion matrix of results for Trial 1 - Constrained 5-stage multigrid partitioning. Labels in the 
rows are the GT labels, and labels in the coiumns are the partitioned labels. 



Labels 


LBL-Ot 


LBL-02 


LBL-03 


LBL-04 


LBL-05 


LBL-06 
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0 
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0 
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6 
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Table X. Confusion matrix of results for TrM 8 ^ Constrained 5 stage multigrid partitioning using a compact 
5x5 neighborhood. Labels in the rows are the GT labels^ and labels in the colunins are the partitioned labels. 



Labels 
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0 


0 


0 


G4 


0 


0 


0 
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u 


0 


0 


G5 


. 0 


0 


0 


0 


172 


0 


0 


G6 


0 


0 


0 


40 


146 


0 


0 


Total GT 


2070 


1602 


2773 


2332 


564 


113 


10 


Total non-GT 


25251 


29977 


47107 


26408 


17295 


17372 


7126 
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Table XL Confusion matrix of results for Trial 9 r Constrained 6 stage multigrid partitioning using a 
compact 5x5 neighborhood. Labels in the rows are the GT labels, and labels in the columns are the partitioned 
labels. 
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E3 
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0 
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o 


A 


0 
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E4 


0 
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15 


0 


0 


0 
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U 
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4" 


A 


A 
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1 
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1 


0 
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J2 


ft 
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U 
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u 


p3 
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0 


1 

1 


A 


Nl 
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0 
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R3 


50 


A 
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U 
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U 


A 


0 


0 


Dl 
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1 
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ft 
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V 


' ft 


1 

1 
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D3 
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7 
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5 


0 


0 


D4 


4 


8 


28 


18 
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18 


0 


D5 


5 


12 


396 


9 


0 


8 


0 


D6a 


10 


• 3 


514 


5 


17 


10 


0 


D6b 


12 


13 


600 


12 


6 


10 


0 


D6c 


2 


0 
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0 


0 


0 


0 


Gl 


0 


1 
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0 


0 


0 


2 


G2 


0 


1055 


0 


1 


' 0 


0 


0 


G3 


2 


0 


0 


0 


0 
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0 


HI 


0 


0 


0 


33 


0 


0 


0 


Wl 


0 


0 


0 


66 


0 


0 


0 


G4 


0 


12 


0 


0 
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2 


b 


G5 


0 


172 


0 


0 


0 


0 


0 


G6 


0 


0 


0 


147 


39 


0 


0 


Total GT 


529 


1718 


2726 


678 


2542 


1261 


10 


Total non-GT 


18699 


36006 


32998 


25531 


29256 


20920 


7126 
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Table XIL Errors of commission for eight trials Each column lists the percentages of labels assigned to a 
partition belonging to a GT category other than the one listed at the top of a column. 



Commission Errors 



Trial 


Trees 


Grass 


' Asphalt 


ISOPATA 


22.48% 


1.32% 


68.66% 


Trial 1 


4.90% 


2.16% 


33.07% 


Trial 2 


• 24.08% 


3.88% 


20,40% 


Trial 5 


11.81% 


12.11% 


86.27% 


Trial 6 


6.09% 


3.43% 


20.19% 


Trial? . 


0.94% 


3.03% 


21.12% 


Trial 8 


21.45% 


34.08% 


81.32% 


Trial 9 


39.03% 


27.82% 


80.26% 
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Table XIIL Measures of clutter, l^ciutter listed for eight of the trials. 



Trial 


^Clutter 


ISODATA 


0.2263 


Trial 1 


0.0170 


Trial 2 


0.0320 


Trials 


0.0720 


Trial 6 


0.0334 


Trial? 


0.0424 


Trials 


0.0401 


Trial 9 


0.0483 



Table XIV. Measures of speckle. Two smooth measures of speckle, K^^^^^^^ and ^%eckie ' ^'^^^^ 
the trials, along with the ratio of the two measures. 



Trial 


K'(/i=2) 
^Speckle 


|^(«=4) 
^Speckle 


^Speckle ' ^Speckle 


ISODATA 


1641 


3359 


2.047 


Trial 1 


222 


479 


2.158 


Trial 2 


395 


795 


2.013 


Trial 5 


665 


1377 


2.071 


Trial 6 


372 


794 


2.134 


Trial 7 


536 


1093 


2.039 


Trials 


467 


930 


1.991 


Trial 9 


477 


924 


1.937 
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Figure 1 A neighborhood system Q^stt site s. 



Figures 




Figures. ISODATA duster map. 




Figure 6. Trial 3 - Unconstrained method using classical KS measure 




Figure 9. Trial 6 ~ Unconstrained method initialized by ISODATA using spectral angle and Eiiclidean 
distance measures. 



Figur^ 10, Trial 7 - Constrained 5-stage method using two disparity measures. 





Figure 12. Trial 9 - Constrained 6-stage method using two disparity measures but a compact 
neighborhood. 
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iSpectral/spatial annealing of hyperspectfal imagery initialized by a 

supervised classification method 

Robert S. Rand 

U.S. Amay Engineer Research and Development Center, Alexandria, VA 223 1 5-3864 



ABSTRACT 

A simulated annealing method of partitioning hyperspectral imagery, initialized by a supervised classification method, is 
investigated to provide spatially smooth class labeling for terrain mapping applications. The method is used to obtain an 
estimate of the mode a Gibbs distribution defined over a symmetric spatial neighborhood system that is based on an 
energy function characterizing spectral disparities in Euclidean distance and spectral angle. Experiments are conducted 
on a 2 1 0-band HYpICE scene that contains a diverse range of terrain features and that is supported with ground truth. 
Both visual and quantitative results demonstrate a clear benefit of this method as compared to spectral-pnly supervised 
classification or unsupervised annealing that has been initialized randomly. 

Keywords: Gibbs distribution, Bayesian, Markov Random Field, supervised classification, Euclidean distance, Spectral 
Angle, clustering, hyperspectral imagery. 

1. INTRODUCTION 

A Bayesian spectral/spatial approach to segmenting hyperspectral imagery into homogeneous regions is investigated 
where spatial consistency is imposed on the spectral content of sites in the scene. The intention is to produce a spatially- 
smooth labeling of regions in a scene with minimal loss of significant information and scene content. For example, we 
desire to eliminate as much insignificant clutter as possible (such as the mottled tree patterns in forested regions) while 
retaining the small, but significant, terrain features (i.e. a small rock outcrop, or building, etc.). 

12 

In prior work, * we investigated a Bayesian Gibbs-based unsupervised approach to partitioning hyperspectral imagery 
■ that is initialized randomly. Various energy ftinctions using different disparity measures were investigated, the merit of 
using a non traditional extended neighborhood system and a Renormalization Group (RB) based multigrid 
implementation to improve global labeling was demonstrated, and an enhanced method of spectral mixture analysis 
(SMA) making use of such a partitioning process was proposed that is appropriate for terrain mapping when the imagery 
involved is characterized by moderate-to-high scene complexity. Whereas, the results have been very encouraging as 

compared to other existing unsupervised methods, such as KMEANS and I SO DATA, '^'^ concerns remain about 
obtaining a truly global labeling solution for practical computing times, and there has been a nagging problem with the 
method having a hard time extracting sniall narrow terrain features, such as roads. However, as we have repeatedly 
pointed out, we only examined its behavior for the case of random initialization, which is really a *'worst-case" scenario. 

The current effort extends the development of the unsupervised spectral/spatial approach by initializing a partitioning 
algorithm with another Bayesian-based classification process. By changing the assumptions governing the class- 
conditional and prior probability distributions we recover a traditional spectral -only supervised process that can provide 
spectrally-optimal initial estimates for the spectral/spatial partitioning process. Initializing the partitioning algorithm in 
this manner effectively makes the method a supervised classifier. 

The prior and current work has its foundations based on the pioneering efforts of others in the Bayesian community. A 
general mathematical development of the underlying theory, applying the Bayesian framework in the context of image 

processing, is well established.^ Using such a framework, the properties of a Markov Random Field (MRF) have been 



applied to develop a Gibbs-based partitioning algorithm for gray-scale imagery with two spatial coordinates, and the 
Gibbs sampler and Metropolis methods have been used as ways to compute maximum a posteriori (MAP) estimates. ' ' 

2. DESCRIPTION OF ALGORITHM 

2.1 Bay esian Framework 

A Bayesian approach provides our framework to construct an algorithm for the classification labeling of homogeneous 
regions in hyperspectral imagery. Consider the combined process^ = where is a discrete labeling 

process that provides class labels that identify similar regions in a scene, and is a spectral modeling process for the 
underlying spectra of materials in a scene of hyperspectraf data. What is actually observed is G, the observed 
hyperspectral image cube. For X^ and <? there are both spatial and spectral coordinates. 

The algorithm is geiieralized to be multigrid, and thus , we define a collection of lattices that provides a finite 
sequence of coarser and coarser grids: ^p'''^ ^ 5^,^'^ ^ ... 

S';''^{{ia,^\J(T„+\): 0<i<%-i,0^y5^} (l) . 

N^^ is the number of stages in the sequence, a„ is a parameter called the "grid resolution" which determines the spatial 
sampling of the algorithm and the spatial resolution (coarseness) of the labeling process at stage w (n = \y..^yN^ ) of the 

multigrid processing, iV^„^^ is the number of rows in an image, and A^^o/s number of columns. S^p"^ is the set of 
sites that will be used in the partition process at a resolution a„. Our grid structure is a special case of the general 
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course-tp-fme grid structure used in the RGA approach to image processing, where the coarsest level contains all the 
large-scale features and each finer level contains correspondingly smaller-scale features. The term "site" refers to a 
generic element of a lattice, as specified in (I). At full-resolution a site corresponds to a pixel, but at a coarser scale it 
corresponds to a square block of pixels. The value of a„- is defined to be the number of pixels per side of each block of 
pixels. The finest resolution is <t, =1 . We denote the full-resoliition image pixel lattice as S,, i.e., the resolution at 
which the observations are obtained, where =Sp^. 

From a Bayesian perspective, the classification problem can be viewed as maximizing the posterior probability 
distribution I ^''^ ), which can be written 

■ ?t{x'' \X^)oc.Pr{x^ \ X')?r{x'') (2) 

where ?r(^X^ \^^) is the class conditional distribution and Pr(^X^ ) is the prior distribution. X^ is a discrete labeling 
process which associates a label with each site. The realizations of are e F, where r = {l,2,..., A^^^^/J is the 
set of classification labels. N^^^^f^ is the total number of labels assigned to the process. The configuration space for this 
labeling process X^ is then determined , by N^^^j^ and a„ (the label resolution), and it is denoted by 

Q<r-' = {x" ^(xf, ^ e S^'') I x^e {h.-.N^,}} • 



2.2 The Partitioning Algorithm 



We now make the assumptipn of iincorrupted data such that^^ =G and of singular correspondence between data 
assignments and class membership P((7 = g | A''' = x'']|= 1 , where g are the realizations of G and x'^are the realizations 

of . Therefore, (2) simplifies to ?v{x^ I ^^)*^ Pr(A"^)i and the classification problem reduces to maximizing 

Pr ( A^'* ) i Under these assumptions we refer to the labeling process as the "Partition Process." 

Furthermore, we treat A''' as a Markov Random Field (MRF) with respect to a graph {sf"\Q^^''^^ , where is the 

entire neighborhood structure for the partition process at resolution g„ , and (g^ is the specific neighborhood at site s for 
the partition process at resolution <T„. For any S that is one of S^p"^ : 5 = {5i,...,^^ }, where Af is the number of sites on 
the grid, an associated neighborhood systeim ^^*'"V={g|j,*e 5} is specified which is a set of lattice graph structures of 
neighbors, such that s^Q^, and se^Q^ if and only if r e . 

In constructing an algorithm to maximize Pr(A'''), the development is similar to that of Geman and Gemah^ with the 

important additions of a spectral dimension and the multigrid generalization. We exploit the Gibbs equivalence 

7 • ' 

theorem which recognizes the MRF as a random process with a Gibbs distribution. The advantage of adopting this 

approach is that we can now compute the maximum a posteriori (MAP) estimate that maximizes the probability 
Pr(Ar^) by iteratively sampling from the local Gibbs distribution pertaining to each site 5g S^p"^ . Our algorithm, called 
the "partitioning algorithm,'* seeks to maximize the distribution: ,y 

• • (3) • 

■ 

where {x^ = x^} means {x^^ = jc;J^, , _X^^ ~ jc^ } , 2, = ^-7^*** ^ ^nd ^/^(jc^,g) is the energy interaction of 

- ' /' " x^kr . * • 

site Sf"^ with the neighborhood Q^. The realizations are estimates of x^ obtained by using the results of the 
previous iteration. Quite conveniently, the computations of Z, are not necessary for purposes of obtaining a MAP 
estimate. We denote g to be the spatial and spectral realizations of G, which is in contrast to the spatial-only realizations 
of the simpler G considered by Geman and Geman.^ Specifically, are the spectral observations at site and g are the 
collections of these observations for all the sites. 

The energy function U^(x^,g) depends on difference measures, called "disparity measures," which indicate the 

dissimilarity between a site 5 and a neighbor r. Our specific energy function combines two disparity measures that 

measure spectral difference in two very different ways: the spectral angle and the Euclidean distance.* This energy 
function is designed to respond to disparities as measured by the spectral angle measure and the Euclidean distance. Our 
reason for using a combined measure is that we seek to increase the sensitivity of the algorithm without introducing 
■undesirable clutter. Rather than using a single measure with a very sensitive threshold that would respond to one 
specific type of spectral difference, we use two measures that respond to different kinds of spectral difference with less 
sensitivity. 

A simulated annealing technique is used to obtain the MAP estimate via the Metropolis algorithm, with an 

eventual fast freeze at low temperature.. The computations involve a temperature at iteration k that defines an 



annealing schedule. Furthemaore, the multigrid implementation uses a grid-coarsening and cascading technique that is 
consistent with aspects of the Renormalization Group Algorithm (RGA) method developed for the restoration of 
8 

degraded images. Specifically, the fine-to -coarse sequence, described above, with partition lattices defined by (1) is 
reversed to a coarse-to-fine sequence: jSp""^^^^ ->...—> 5^^'^ ^ S^p^ , The results of stage "w" are used to initialize the 
aigorithm at the next finer stage "w-/." ' 

2.3 Initializing the Algorithm 

The algorithm requires an initial estimate of the realizations. These values can be initialized randomly, in which 
case during the first iteration the estimated realizations are assigned class labels randomly. This is a worst-case 
scenario. Altematively, we can consider some other means of initialization. 

One alternative is to initialize the algorithm with the results of some prior supervised classificafipn run. Consider (2) 
and recall our assumptions X^ = G and P[G = g\X^ = x^)=\ that simplifies (2) to Pr{X'' |jr'^)oc Pr(A^''). This 

was actually a gross simplificafion required to keep the algorithm analytically tractable and computationally practical so 
that we can incorporate spatial information into (2). . 

A common spectral-only Bayesian approach is to use (2) and assume that the siterSpecific distributions are independent 
so that (2) becomes Pr(x^ | Jr^) = P[Pr{A'^| )Pr(vVf ) and to further assume equal prior probabilities (all the 

Pr^X^'j are equal). Consequently, each site iri an image can be classified independently of the other sites by selecting 
jcf which maximizes ?v{^X^ \ X^Y Furthermore, these class-conditional distributions are often assurned to be 
multivairiate normal, which leads to the classification rule: 

if =argmaxVv(jU^,i:^) (4) 

where ju^ is the mean, and is the covariance, of the spectral vectors for the class label x^ , respectively. Eq. (4) is a 

■ 4 5 

well-known quadratic multivariate classifier used in many pattern recognition scenarios. ' Although the multivariate 

normal assumption is often violated, it provides a far better approximation to the true distribution the class-conditional 
spectral data than the simplifying assumpfipn P{.G- g\ X^ =^ x^)=\ made above. 

The classification rule in (4) can be simplified further if we make additional assumptions on the class covariance 
matrices. Specifically, if we make the approximation that = I, the Identity Matrix, for alljcf e F, the classificafion 

4 5 

rule simply reduces to the Euclidean distance classifier, which is a linear classifier. ' 

The partitioning algorithm has the advantage of incorporating spatial information and smoothing but the disadvantage of 
using a crude class-conditional distribution. The classifier provided in (4) uses a better approximation to the true classr 
conditional (spectral-only) distribution and can be shown to be the optimal classifyier when the underlying assumptions 
are met; however, it does not model any spatial interactioti or provide any spatial smoothing that could be provided by a 
spectral/spatial prior distribution Pr(X'*) . 

Initializing the partitioning algorithm with the results of a classifier that computes (4) offers the potential for improved 
classification by providing initial estimates based on label estimates that under the proper conditions are spectrally 
optimal. Accordingly, we propose the following two-step procedure; 



Step 1 : Perform a classification of the scene using a supervised classifier of the form (4). 

Step 2: Perform a spectral/spatial partitioning of the scene based on (3) initialized by the results in Step I. 

Implemented in this manner, the algorithm can also be viewed as a spectrally-optimal supervised classification algorithm 
that has been smoothed by a post-processing routine that . imposes spectral/spatial constraints defined by the Gibbs prior 
probability distribution Pr (>V ) . 

3. DESCRIPTION OF EXPERIMENT 

II. Data \ . 

A HYDICE scene is used for validation that was extracted from data collected over a geographic area near Fort Hood, 
Texas at a 10,000 feet altitude, on 16 May 1996. Figure 1 shows the scene, consisting of 300 samples by 600 lines of 
pixel vector data with a spatial resolution of approximately 1.5 meters. Because of the severe degradations that occur in 
atmospheric windows of the spectruni and sensor degradations occurring in certain bands, only 1 17 out of a possible 210 
band were used to compute the disparities within the energy functions. These 1 17 bands were calibrated to reflectance 
by the Empirical Line Method^ based on the known spectra of the calibratioii panels that were located iti a region outside 
this particular scene. 




Figure 1. HYDICE Scene (Band 49) 



Ground truth, describing certain ground features in the scene, was based on data collected during a site visit to the area 
during April 2000. In addition to the grpund-featiire data, spectral data of a diverse range of materials, measured by field 
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spectrometers and collected independently of this experiment at other geographic locations, were used in prior work * 
to determine an appropriate quantization threshold for the spectral angle and Euclidean distance disparity metrics, which 
remained the same for this v/oik. Ground-measured spectral data within the study site were not available. 

i.Z Validation of the Algorithm ' * 

An experiinent was conducted using computer code written in C++ that iniplements the algorithm described in Section 2. 
Two trials were conducted, investigatiiig the merit of using the two-step supervised approach as compared to simply 
using the randomly initialized annealing approach. 

In Trial A, the randomly initialized partitioning algorithm was run using the parameters shown in Table 1 for the various 
stages (grid sizes) of the multigrid process. We set Ni^i^, =10. As indicated in Table 1', the trials proceeded in a 
multigrid sequence with grid resolutions a = 16, 8, 4, 2, 1, 1 Simulated annealing was used to compute the MAP 



estimates using an energy fiinction comprised of the combined disparity measures that is identical to that used in our 

12 ' ' 

previous work. ' A range C ^ (0.05-3.0) was used, which is below the upper bound, and influences algorithm 

convergence as discussed in this previous work. The value of C was set to its highest value at the initial stage of the run 

(corresponding to high temperature) and it was set to lower values for the subsequent stages. In addition to the practical 

computational advantage, a lower C value allows the algorithm to "remember" the results at previous stages. The 

distances of the intermediate and far neighbors shown in Table 1 corresponds to fixed multiples of the grid resolution, 

4a and 8a, respectively. The algorithm was initialized randomly, and then proceeded from a grid resolution of a^l6 to 

a=l, using the extended neighborhood of 24 neighbors. A final pass was then made at CJ=7 with a compact 

neighborhood comprised of 8 near neighbors. 



Table 1. Control Parameters for Trial A. The sequence of numbers \n the "Dist of neighbors" column 
corresponds to the distance from the "center pixel (the "site") to the perimeter, for the surrounding near 
neighbors, the intermediate neighbors, and the far neighbors. The sequence of numbers in the **# of neighbors" 
column corresponds to the number of near neighbors, intermediate neighbors, and far neighbors in the 
neighborhood system. 



a 


C 


# iterations 


Dist. of neighbors 


# of neighbors 


16 




20,000 


16,64,128 ' 


8,8,8 ' " 


8 


2.0 


2,000 " 


8,32,64 


8, 8, 8 


4 


2.0 


2,000 


4,16;32 


8,8,8 " 


2 


0.5 


350 


2,8,16 


8, 8, 8 


1 


0.5 


350 


1,4,8 


8,8,8 . 


1 


0.05 


300 


1 
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In Trial B, the two-step procedure was implemented. For Step 1, a supervised classification was performed using the 
Euclidean distance classifier. Recall that this classifier satisfies (4), given the assumption that = I for xf g T . We 

use this classifier rather than the quadratic classifier for two reasons: (1) We want to establish the degree of improvement 
over the randomly-initialized partitioning algorithm when using the simplest of assumptions on the class-conditional 
distributions, as well as to ascertain the capability of the partitioning algorithm to "clean up" the results of a rather 
simple classifier; (2) Using a quadratic classifier on the 117-band (very high-dimensional) HYDICE imagery requires 
the intermediate step of applying a dimension reduction technique (such as the Principal Component or the Minimal 
Noise Fraction transformation) to Jhe data which implicitly adds another step that we have not discussed into the 
algorithm. . . 

Using the information obvious in the scene supplemented by ground truth information 14 training classes were extracted 
fi-om the scene. Some of these classes represent similar terrain features, such as Light Asphalt, Medium-toned Asphalt, 
and Dark Asphalt, which were subsequently combined after the classification run was finished. They were combined 
afi:erwards because combining such features into a single training class would either result in a multimodal class- 
conditional distribution or a unimodal class-conditional distribution of extremely high variance. So after the 
classification run was finished, the 14 classes were merged into 9 classes. The merged classes along with the 
corresponding original training classes are listed in Table 2. Further information about these classes is documented in 
ground-truth notes, but not discussed here. 

For Step 2, the partitioning algorithm was tested two different ways. First, the algorithm was initialized by the results of 
Step 1 at a grid size of a=4, and then set to run through the multigrid sequence down to<T=J according to the parameters 
in Table 1. Second, the algorithm was simply initialized at the finest grid size cmI and then run at low temperature 
(C==0.05) with the parameters shown in the bottom row of Table 1 . 



Table 2. Description of Classes. 



Class Name 


Description 


# of samples 


Asphalt 


Light Asphalt, Med. Asphalt, Daik Asphalt 


312 


Concrete 


Bright concrete 


120 


RoofL 


Two light rooftops - Roof Bl, RoofB2 


162 


RoofH 


Rooftop of asphalt-singled houses 


17 


Roof B 


Two brown rooftops - Roof J3, J2 


80 


Soil 


Exposed soil 


56 


TreesDrk 


Dark Green Trees ' 


430 


Grass 


Tv/o grass regions of differing vigor 


727 


GrassStr 


Highly stress grass region 


38 



4. RESULTS 

Figure 2 shows the resulting class map for the randomly initialized unsupervised partitioniiig algorithm in Trial A using 
the parameters in Table 1 . The results are very similar to the results for A^^^^^/, =6 obtained in our previous study. ^ 
Visually, the partitioning algorithm generates smooth partitions that accurately represent structured regions in the image. 
For most of the bigger regions, pixels with the same label represent the same type of phenomenon on the groimd, 
globally across the im^ge, and regions with different labels correspond to different phenomenon. However, a globally 
consistent labeling of the partitions was not completely achieved. In particular, labeling problems can be seen in some 
of the smaller regions. Although the contiguous pixels in smaller partitions represent the same phenomenon, locally, 
some small partitions, separated by a certain amount of distance and given the same names, actually represent different 
phenomenon; Conversely, some small partitions, separated by a certain amount of distance and giyen different names, 
actually represent similar phenonienon. 

Figures 3 and 4 show the resulting class maps for Trial B. The class map in Figure 3 is the result of the spectral-only 
Euclidean classifier. Fourteen training classes were merged into the 9 thematic classes listed in Table 2. Visually, the 
map accurately portrays many of the structures in the scene and does a better job at extracting linear terrain features, 
such as roads, than the randomly initialized algorithm in Trial A; however, it does so with a severe amount of speckling. 
Figure 4 shows the results of the latter of the two implementations tested for Step 2 the implementation for which the 
algorithm was initialized at the finest grid size a^I and then run with the parameters shown in the bottom row of Table 
1. Visually, this class map shows a significant improvement over those shown in Figures 2 and 3. The method retains 
the structures and the linear terrain features extracted by the Euclidean classifier while eliminating most of the speckling. 

Regarding the other niethod tested for Step 2 (results not shown), this implementation tended to lose some of the linear 
features identified by the Euclidean classifier. Although the degradation was not severe, this approach is also 
computationally more expensive. So, whereas a multigrid implementation is the preferred approach when the algorithm 
is initialized randomly, the latter (single pass) implementation is preferred for the two-step method. 

Tables 3 and 4 provide siipporting quantitative results regarding the labeling in the Trial B class maps. Tables 3a and 3b 
provide auto-classificatioti results, listing the percentage of sites in each of the training regions that have been assigned 
to a particular training class, for Step 1 and Step 2, respectively. Tables 4a and 4b provide classification results for test 
data extracted from the scene based on ground-truth, listing the percentage of sites in each of the test regions that have 
been assigned to a particular training class. 



Figure 2. Trial A — Randomly-Initialized Multigrid Partitioning Algorithm. 




Figure 4. Trial B (Step 2) — Partitioning Algorithm Initialized by Supervised Classification. 



Table 3. Quantitative Results for Test Data in Trial B. The columns list the percentage of sites in each of the 
training regions that have been assigned to a particular training class (i.e. Training data are tested against itself). 
Percentages highlighted in bold (and expected to be on the diagonal) are considered to be correct assignments. 

Table 3a. Results for the Euclidian Distance Classifier. 





Asphalt 


Concrete 


RoofL 


RoofH 


RoofB 


Soil 


Trees 


Grass 


GrassStr 


tfptS 


Asphalt 


99.7% 


0.0% 


0.0% 


0.0% 


0.3% 


0.0% 


0.0% 


0.0% 


0.0% 


312 


Concrete 


6.7% 


. 54.2% 


31.7% 


0.0% 


0.0% 


7.5% 


0.0% 


0.0% 


0.0% 


l2u 


RoofL 


1.2% 


4.9% 


92.6% 


1.2% 


0.0% 


0.0% 


0.0%. 


0.0% 


0.0% 


162 


RoofH 


0.0% 


0.0% 


5.9% 


94.1% 


p.0% 


0.0% 


0.0% 


0.0% 


0.0% 


i 7 


RoofB 


6.3% 


0,0% 


0.0% 


0.0% 


93.8% 


0.0% 


0.0% 


0.0% 


0.0% 


■ O/l 

OU 


Soil 


0.0% 


1.8% 


0.0% 


0.0% 


0.0% 


98.2% 


. 0.0% 


0.0% 


0.0% 


56 


TreesDrk 


3.5% 


0.0% 


0.0% 


0.0%. 


5.8% 


0.0% 


63.5% 


27.2% 


0.0% 


4J0 


Grass 


0.1% 


0.0% 


0.0% 


0.0%. 


0.0% 


0.0% 


3.6%. 


95.9% 


0.4% 


/27 


GrassStr 


b.0% 


0.0% 


. 0.0% 


0,0% 


0.0% 


0.0% 


0.0% , 


0.0%* 


100.0% 
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Table 3b. Results for the Partitioning Algorithm Initialized by the Euclidian Distance Classifier. 






Asphalt 


Concrete 


RoofL 


RoofH 


RoofB 


Soil 


Trees 


Grass 


GrassStr 


Upts 


Asphalt 


100.0% 


0.0% 


0.0%' 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


312 


Concrete 


0.8% 


86.7% 


12.5% 


0.0% 


O.OVo 


0.0% 


0.0% 


0.0%, 


0,0% 


120 


RoofL 


0.0% 


. 0.0% 


100.0% 


0.0% 


0.0% 


Q.0% 


0.0% 


0.0% 


0.0% 


162 


RoofH 


. 0.0% 


0.0% 


0.0% 


100.0% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


17 


RoofB 


0.0% 


0.0% 


0.0% 


0.0% 


100.0% 


0.0% 


0.0% 


0.0% 


0.0%- 


80 


Soil 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


100.0% 


0.0% 


0,0% 


0.0% 


56 


TreesDrk 


2.1% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


82.8% 


15.1% 


0.0% 


430 


iGrass 


0.0% 


0-0% 


0.0% 


0.0% 


,0.0% 


0.0% 


0.7% 


993% 


0.0% 


727 


GrassStr 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


100.0% 
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These tables demonstrate a clear benefit to using the algorithm. Table 3 shows that when the algorithm is applied to its 
own training data, the percentage of correct labeling increases for eight out of the nine classes. For the ninth class, the 
classification accuracy was already 100% so there was no possibility of improvement. Table 4 shows that when the 
algorithm is applied to test data from regions outside of the training data there is cotisistent and significant improvement 
in accuracy. Out of the 15 test regions, there was improvement in nine of the regions, there was no improvement for 
three of the regions, there was ivp chance of improvement for two of the region (100% accuracy), and a decrease in one 
of the regions. Notice that the degree of improvement was sometimes substantial. For example, the accuracy of the test 
region for trees improved from 66.9% to 87.9%. 

In the case where there was a decrease in accuracy, the partitioning algorithm was still behaving in an appropriate 
manner. The spectral signatures of concrete and soil are quite similar (concrete being composed of soil/silicate 
materials) and consequently simple classifiers often confuse these two classes. In fairness to the algorithm, in this case, 
the Euclidean classifier incorrectly labeled 66.7% of the soil region as concrete and another 10.7% as asphalt, leaving 
only 20.1% of the sites correctly labeled for the initialization step. Understandably, the partitioning algorithm then 
proceeded to expand the incorrect concrete labeling to 91.8%. 



Table 4. Quantitative Results for Test Data in Trial B. The columns Ust the percentage of sites in each of the test 
regions that have been assigned to a particular training class. Percentages highlighted in bold are considered to be 
correct assignments. ' . . 



Table 4a. Results for the Euclidian Distance Classifier. 





Asphalt 


Concrete 


Roof L 


RoofH 


Roof B 


Soil 


Trees 


Grass 


GrassStr 


Upts 


AsphaltLgt_T 


80,4% 


0.0% 


0.0% 


0.0% 


18.5% 


0.0% 


0.0% 


Ll% 


0.0% 


92 


AsphaltMed_T 


86.6% 


0.0% 


0.0% 


0.0% 


13.4% 


0.0% 


,0.0% 


0,0% 


0.0% 


.97 


AsphaltDrkT 


100.0% 


0.0% 


0.0% 


0.0% 


. 0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


437 


Rooffll T 


0.0% 


0.0% 


100.0% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


. 0.0% 


173 


Goncrete_T 


0.0% 


87.5% 


12.5% 


0.0% 


0.0% 


0.0% 


0;0% 


0.0% 


o:o% 


24 


RoofGravel 


^ 11.8% 


82.4% 


5.9% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


17 


RoofB2_T 


6.3% 


17.2% 


76.6% 


0.0% 


0.0% 


0.0% 


0.0% 


; 0.0% 


0.0% 


fid 


TennisCburt 


2.3% 


0.0% 


- 0.0% 


97.7% 


0.0% 


0.0% 


0.0% - 


. 0.0% 


' 0.6% 




RoofHpuse_T 


0.0% 


0,0% 


24.4% 


66.7% 


■ 8.9% 


0,0% 


0.0% 


0.0% 


0.0% 




RooD3^T 


5.1% 


0.0% 


0.0% 


0.0% 


94,9% 


0.0% 


0.0% 


0,0% 


0.0% 




Roofl2_T 


7.1% 


0.0% 


0.0% 


0.0% 


92.9% 


0.0% 


0,0% 


0.0% 


.0.0% 


A") 


SoilLgtT 


0,0% 


11.6% 


^ 0.0% 


0,0% 


0.0% 


88.4% 


0.0% 


0.0% 


. 0.0% 


1 JO 


SoilTan_T 


10.7% 


66.7% 


0.6% 


0.0% 


0.0% 


20.1% 


0.0% 


0.6% 


1.3% 


1 Jy 


GrassWild^T 


0.2% 


0,0% 


0.0% 


0.0% 


0.7% 


0.0%. 


18.9%' 


74,7% 


.5.5% 


Jinn 


TreesDrk_T 


5.5% 


0.0% 


0.0% 


0,0% 


11.6% 


0.0% 


66.9% 


16.0% 


0.0% 


jjyj 


GrassStr_T 


0.0% 


0.0% 


6.0% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


100.0% 




Table 4b. Results for the Partitioning Algorithm Initialized by the Euclidian Distance Classifier. 






Asphalt 


Concrete 


RoofL 


RoofH 


RoofB 


Soil 


trees 


Grass 


GrassStr 




AsphaltLgt_T 


84.8% 


0.0% 


0.0% 


0.0% 


14.1% 


0.0% 


0.0% 


i.1% 


0.0% 


92 


AsphaItMed_T 


96.9% 


0.0% 


0.0% 


0.0% 


3.1% 


0.0% 


0.0% 


0.0% 


0.0% 


97 


AsphaltDrkT 


100,0% 


0.0% 


o:o% 


0.0% 


0.0% 


0,0% 


0,0% 


0.0% 


0.0% 


437 


RoofBl T 


0.0% 


. 0.0% 


100.0% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


173 


Concrete_T 


8.3% 


87.5% 


4.2% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


24 


RoofGravel 


0.0% 


100.0% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


17 


RoofB2_T 


0.0% 


. 0.0% 


100.0% 


0.0% 


0,0% 


0.0% 


0.0% 


0,0% 


0.0% 


64 


TennisCourt 


0.0% 


0.0% 


0.0% 


100.0% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


43 


RoofHouse_T 


0.0% 


. 0.0% 


22.2% 


77.8% 


0.0% 


0.0% 


0.0% 


0.0% 


0.0% 


45 


RoofJ3_T 


6.0% 


0.0% 


0.0% 


0,0% 


100.0% 


0.0% 


0.0% 


0.0% 


0.0% 


39 


Roofr2_T 


0.0% 


. 0.0% 


0.0% 


'0.0% 


100.0% 


0.0% 


0.0% 


0.0% 


0.0% 


42 


SoilLgtT 


0.0% 


n.6% 


0.0% 


0.0% 


0.0% 


88.4% 


0.0% 


0.0% 


0.0% 


138 


SoilTan_T 


3.8% 


91,8% 


0.0% 


0.0% 


0.0% 


. 1.9% 


0.0% 


0.6% 


1.9% 


159 


GrassWiid_T 


0.0% 


0.0% 


0.0% 


0.0% 


0,0% 


0.0% 


20.2% 


74.8% 


5,0% 


1700 


TreesDrk_T 


3.0% 


0.0% 


0.0% 


0.0% 


3.i% 


0.0% 


87.9% 


6.0% 


0.0% 


1591 


GrassStrjr 


0.0% 


0.0% 


0.0% 


0.0% 


.0.0% 


0.0% 


0.0% 


0.0% 


100.0% 


35 



The results show that spatially smooth labeling can be achieved without decreasing the accuracy of classification. 
Tables 3 and 4 show an overall iiicrease, not decrease, in classification accuracy when going from the spectral-only to 
the spectral/spatial process. Furthermore, Figure 4 shows that much of the speckling and edge artifacts in a scene 
labeling can be removed without removing spectrally si^ificant individual pixels. For example, observing the parking 
lot adjacent to the department store (upper-left corner) in Figure 4, the single vehicles are not removed, whereas the 
edge-artifacts of the roads and the clutter in the forbst and grass regions throughput the scene are mostly removed. 

S. CONCLUSIONS 

A Bayesian framework is used to develop a 2-step supervised classification algorithm that is capable of perforrriing high- 
quality spatially smooth labeling of hyperspectral imagery. To demonstrate the concept, a linear classifier was used to 
initialize a Gibbs-based partitioning algorithm resulting in significantly improved label accuracy and smoothness as 
compared to the linear classifier itself. In addition, the global labeling accuracy was increased as compared to the stand- 
alone randomly initialized partitioning algorithm. Based on optimality arguments, the current results will likely improve 
further, if the linear classifier is replaced with a quadratic classifier. 

In order to achieve a reasonable global labeling accuracy, the randomly initialized (unsupervised) partitioning algorithm 
is best implemented as a multigrid process using an extended neighborhood systenri. However, initializing the 
partitioning algorithm with even the simplest of spectral-only supervised classifiers not oi>ly improves the global 
accuracy, but it also reduces the computation by (1) eliminating the need for a multigrid process arid (2) allowing the 

annealing to start at a cooler temperature. 

The experiment shows that spatially smooth labeling does not have to occur at the cost of classification accuracy, which 
has been the case with a number of post-processiiig nriethods that simply app]y spatial constraints to a spectral-based 
classification. ' - ■ 
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